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ABSTRACT 

The solution of a beam on an elastic nonlinear foundation by the 
Finite Element Method is presented in this thesis. Discontinuous 
(Winkler) and continuous foundations are considered. The general 
formulation is based on the Galerkin method. The solution technique 
for the linear case uses Gauss elimination and the solution technique 
for the nonlinear case is based on Brown*s method (a modified Newton- 
Raphson). Some illustrative examples are presented. A brief com- 
parison of continuous and non- continuous foundations is made. 
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I. INTRODUCTION 



In treating the various problems related to structures on elastic 
foundations, it has been customary for engineers to assume that the 
pressure in the foundation is linearly proportional at every point to the 
deflection of the foundation at that point, independent of the pressure or 
deflection occurring in other parts of the foundation. This assumption, 
is mathematically by far the simplest that one can make regarding the 
nature of a supporting elastic medium. It assumes a complete lack of 
continuity in the material of the foundation as if it consisted of a series 
of independent springs which deflect when directly loaded. This theory 
has proved adequate for calculating stresses and deflections in railroad 
tracks and has found many applications to plate and shell structures 
[Ref. l]. However no particular claim has been made that the deforma- 
tion or pressure distribution in actual earth foundations could be pre- 
dicted by this method. 

In recent years, some authors have treated foundations as a con- 
tinuous medium which, in contrast to the Winkler hypothesis, repre- 
sents the case of complete continuity in the material [Ref. 9]. It is 
evident that the two types of foundation models lead to different results, 
but how quantitatively significant is the difference has not previously 
been made clear. 
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The properties of actual foundations, however, seem to lie some- 
where in the gap between these two extreme cases. The physical prop- 
erties of soils are obviously of a very complicated nature. Their de- 
formation behavior is influenced by a number of factors such as the 
physical structure, porosity, existence and movement of fluids in the 
pores, etc. In addition, such geologic features as faults, joints, seams, 
crushed zones, fissures and other tectomic effects produce behavior 
significantly different from that derived on the assumption of continuous 
mass [Ref. 14]. The question is then raised; does the approximation 
based on the Winkler hypothesis still hold over various types of soil 
foundations? If it does not, then how much error may be expected? 

How much does the shearing effect in the material of the foundation 
contribute to its deformation? Considering only the deflections of one- 
dimensional beams on foundations (for simplicity), this thesis attempts 
to answer these questions. 
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II. DESCRIPTION OF ELASTIC FOUNDATIONS 



This chapter presents brief descriptions of the Winkler, Modified 
Winkler and Continuous foundations considered in this thesis. Many 
other foundations have been proposed over the years by numerous 
investigators. Kerr [Ref. 31] presents a summary of a number of 
linear foundations. 

A. DEFINITION OF WINKLER, MODIFIED WINKLER AND 

CONTINOUS FOUNDATIONS 

The following definitions will be used throughout this thesis. 

1 . Classical Winkler Foundation 

This type of foundation is characterized by the fact that the 
deflection at every point of the foundation is linearly proportional to the 
pressure applied at that point and is independent of pressures or deflec- 
tions acting elsewhere in the foundation. This assumption is equivalent 
to considering the foundation as a discontinuous medium composed of 
independent elastic springs. It is believed that this hypothesis was 
proposed in 1867 by Winkler [Refs. 27, 9, ^9]. 

2. Modified Winkler Foundation 

This type of foundation, still a discontinuous medium, is 

^Some authors claim that Euler seems to have been the first to 
formulate the hypothesis, although it is generally attributed to Winkler. 
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however characterized by the fact that the pressure at every point is 
nonlinearly proportional to the deflection at that point. 

In this thesis, reference to a Winkler foundation means either 
the classical (linear) Winkler foundation or the modified (nonlinear) 
Winkler foundation. 

3. Continuous Foundation 

In contrast to the Winkler foundation, the pressure and deflec- 
tion at a point in a continuous foundation is affected by the behavior of 
the entire foundation [Refs. 1, 2?]. 

4. Mixed Mode Foundation 

The real foundations of natural soils are more complex than 
either the Winkler or continuous foundations. They are neither com- 
pletely continuous nor discontinuous [Refs. 14, 6]. The pressure and 
deflection at any point are nonlinearly dependent. This type of founda- 
tion is therefore considered as a combination of the modified Winkler 
foundation and the continuous foundation. 

B. MATHEMATICAL MODELS 

In this chapter the governing equations of the foundations considered 
in this thesis are developed. 

Consider an elastic foundation which is a compressible single layer 
of thickness H placed on a rigid base. It will be assumed that the thick- 
ness of the elastic foundation, its support conditions, the elastic con- 
straints and all other properties as well as the external load do not 



15 



vary in the z_direction. A narrow plate of thickness w is cut from the 
elastic foundation by two planes parallel to the x y plane. A beam of 
the same thickness w rests on the surface of this elastic foundation. 

lateral loading 




Figure 1. Geometry of an elastically supported beam 



Let an external load q(x), pounds per unit length, act on the beam 
(Figure 1). If r(x) is the reaction of the elastic foundation against the 
beam, the differential equation of bending for the beam is then [Ref. 8]: 
(EIV")" = q(x) - r (x) (1) 
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where El (x) and V (x) are the flexural rigidity and the deflection of 
the beam, respectively. 

Equation (1) is true for any type of foundation. It contains two 
unknown functions, V (x) and r (x). In order to determine them, an 
additional relationship is needed. This relationship is associated with 
the response of the elastic foundation and depends on the type of 
foundation being considered. 

1 . Winkler Foundations 

a. Governing Equation 

The Winkler foundation is characterized by a foundation 
modulus, generally denoted by the letter kt Let V, inches and r, pounds 
per inch, be the deflection and the reaction of the foundation, respectively. 
The Winkler foundation is characterized by the equation 

r (x) = k*vP{x) (2) 

For the classical Winkler foundations, P = 1 and k*is a 
constant and is expressed in pounds per cubic inch. This assumption 
has proved adequate for calculating stresses and deflection in railroad 
tracks [Ref. l]. 

b. Values of Foundation Modulii 

For sand, k*may vary from 25 to 100. Ib/cu. in. For 
ordinary soils on which railroad tracks have been built, k*may vary 
from 110. to 130. Ib/cu.in. For gravel, the value of is even more 
uncertain and it can vary from 200, to 1200. Ib/cu. in. [Refs. 1, 27]. 
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2, Modified Winkler Foundations 



For actual foundations of natural soils, k*is the coefficient 
associated with the nonlinear term [Ref. 6]. Therefore, to generalize 
the Winkler hypothesis, let p be any positive real number. Such 
foundations are called modified Winkler foundations. 

Continuous Foundations 

a. Governing Equation 

Let s be the shear force in the foundation, pounds per 
inch, and l^be the foundation modulus, pounds per cubic inch. The 
shear force s is a measure of the friction force between soil particles. 
The governing equation of continuous foundations is [Ref. 9]: 

r (x) = k*V (x) - (s v'(x))' (3) 

Details of the continuous foundation analysis are given in 

Appendix A. 

b. Values of s and k 

The values of s and k depend on the contact area of beam 
and foundation per unit length of beam. For a loaded area of 6. 5 inches 
of width, equal to that of a 12WF2? standard beam, and 1. inch of length, 
typical values of s and k for some soils are: 

Average sea floor sediment: 
s = 18. 6 lb k = 8. 3 Ib/in^ 

A typical beach soil: 

s = 1000. lb k = 500. Ib/in^ 

where k = wk* 
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r> 






A typical inland soil: 
s = 6000. lb k = 3000. Ib/in^ 

Additional values of soil characteristics are given in 

Reference 7. 

4. Mixed Mode Foundations 
a. Governing Equation 

The governing equation of this type of foundations can be 
derived from the definition of the mixed mode foundations given in 
part A. The mixed mode foundation is a combination of a modified 
Winkler foundation and a continuous foundation: 

r (x) = k vP(x) - (s v'(x))' (4) 

It is seen that any of the previous foundations are special 
cases of the mixed mode foundation. For the Winkler foundation s = 0 
and P = 1, for the nonlinear foundation s = 0 and for the continuous 
foundation P = 1. 
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111 . GENERAL FORMULATION OF THE PROBLEM 



A._ GOVERNING EQUATIONS 

Consider a beam on an elastic foundation. The differential 
equation of bending for a beam on an elastic foundation is given by 
[Ref. 8]: 

(Eiv'*(x)) = q(x) - r(x) (1) 

where q(x) is the external loading, V(x) is the beam deflection and r(x) 
is the reaction of the foundation. The mixed mode foundation, by 
definition, can be considered as the most general type of foundation. 
Its governing equation is given Eqn. 1 with 

r(x) = kvP (x) - (sv'(x))' (4) 

and p > 0. 

Substitution of (4) into (1) gives: 

(Eiv"(x))" = q(x) - kvP (x) + (sv'(x)) (5) 

or 

(EIV) - (sV) + kV^ = q (6) 

To solve this equation for any positive real p , consider a least 
square best fit polynomial such that 

P 2 

V (x) = bj + b^ V(x) + b^ V^(x) + . . . 
where the b^ coefficients (i = 1,2,3,...) are constants. 



( 7 ) 



For simplicity in the development of an actual computer program, 
equation (7) is taken to be a second degree polynomial. This restriction 
is not inherent in the Galerkin procedure but rather in the Finite Element 
formulation of Galerkin adopted in this thesis. Substitution of (7) into 
(6) leads to the following equation: 

(EIV)" - (sV)' + k (bi + b 2 V + b^V^) = q (8) 

or 

(EIV")'' - (sV)' + kb 2 V + kb^V^ = q _ kbj (9) 

Letting 

aj = kb^ (10) 

a2 = kbj (11) 

f = q(x) - kb^ (12) 

equation (9) may be rewritten; 

(EIV")" _ (sV)' + ajV + a 2 V^ = f (13) 

1 . Beam on a Winkler Foundation 

Equation (13) can be applied to the case of a beam on a classical 
Winkler foundation, by letting bj = b^ = 0, b 2 = 1 and Sj = 0. Hence 
the governing equation of bending for a beam on a classical Winkler 
foundation is : 

(EIV)’* + a^V = q(x) (14) 

where a^ is defined by Eqn. 10. 

2. Beam on a Modified Winkler P'oundation 

When s = 0, Equation (13) becomes the governing equation 
of bending for a beam on a modified Winkler foundation: 
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f 



(15) 



(EIV)" + ajV + = 

where terms are defined by Eqns. 10, 11, 12. 

3. Beam on a Continuous Foundation 

Equation (13) can be applied to the case of bending of a beam 
on a continuous foundation, by letting ^2 ” ^1 ” ^3 “ ^ * 

(EIV")" - (sV*)* + ajV = q(x) (16) 

4. Beam on a Mixed Mode Foundation 

Equation (13) is the governing equation of bending for a beam 
on a mixed mode foundation: 

(EIV)" - (sV)' + ajV + = f (17) 

B. BOUNDARY CONDITIONS 

Reference 8 thoroughly develops the boundary conditions for a beam. 
Only the boundary conditions on displacement and slope (the so called 
essential or principal boundary conditions) must be imposed in a 
variational formulation of the problem. 

C. ASSUMPTIONS AND RESTRICTIONS 

Some restrictions have already been stated in the introduction 
section of this thesis; other restrictions will now be discussed. 

1. The flexural rigidity El must be "slowly varying" in accord- 
ance with the development of the basic equations of classical beam 
theory. 
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2 , 



The shear coefficient s must also be "slowly varying" in 



accordance with the development in Appendix A. 

3. The foundation modulus k(x) and the load variable q(x) need not 
be slowly varying. 

D. SOLUTION TECHNIQUES 

The Finite Element Method via the Galerkin approach will be used 
here. Its formulation will be presented in the next section. There are 
two types of problems to be considered. 

1 . Linear problems 
Consider the linear system 

{EIV")'' . (sV»)* + kV = q (18) 

Suppose the solution for q = q^^ is and the solution for q = 

^2 9 ~ ^ ^2 have the solution: 

V = Vj + a (19) 

The application of the Finite Element Method to any linear 
differential equation leads to a linear system of algebraic equations 
which will be solved by Gaussian elimination. 

2. Nonlinear problems 

In contrast to the linear case, any nonlinear differential 
equation, like (13) or (15), leads to a nonlinear system of algebraic 
equations which, in this thesis, will be solved by Brown's iteration 
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method. Reference 13 thoroughly develops the Brown's iteration method 
which is a modified Newton-Raphson method. 

The results of the Finite Element Method solutions presented 
here will be checked with the classical solutions if available, or with 
other methods, such as the Finite Difference Method. Programming 
details will be presented in Chapter V. 
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IV. FINITE ELEMENT FORMULATION 



A. APPROXIMATION BASED ON THE GALERKIN PROCEDURE 

In order to obtain a numericaJ .solution for displacements of beams 
on elastic foundations, let V be approximated by the m degrees of 
freedom expression: 



m 



V “ E '^i <*> 

i= 1 



( 21 ) 



where G^(x) are global (or system) shape functions and the unknowns 
coefficients are generalized coordinates [Ref. 14], i. e. , system 
degrees of freedom. Since the Galerkin method deals directly with the 
differential equation [Refs. 21, 15, 22], (21) will be inserted in (13) 
and the equation residual is formed as follows: 

[<(EIV.G. )”, Gj^(x) > - < (s V.G.'(x)), Gj^(x)> 

i= 1 

m 

+ <aj V.G.(x), Gi,(x) > + ^ <a V. G. (x)V G (x), Gj^(x) > ] 

3=1 1 J J 

- < f, Gj^(x) > = 0 

where the notation < a,b> means jQa(x)b(x)dx . Then 

[V. < (EIG!'(x)) , Gj^(x) > - V. < (s g!(x)), Gj^(x) > 



i=l 



m 



+ a^V. < G.(x), Gj^(x) > + ^ < Gi(x)Gj(x), Gj^(x) > ] 

j=l 

- < f, G. (x) > =0 
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Integration of the first two terms by parts yields: 



V [V, < eig!'(x), gm > + V. < sg!(x), g' (x) > 

1 1 1 1 K 

1=1 

m 

+ a^V. < G.(x), Gj^(x)> + ^ G.(x)Gj(x),Gj^(x) > ] 

- < f, Gj^(x)> =0 (22) 

There is a constant of integration due to boundary conditions left 
by the process of integration by parts in Eqn. 22. However this constant 
can be omitted in the Galerkin procedure [Ref. 23], 



m 



i=l 



Let 










= 


<EIg|'(x), 


Gj^'(x) > 


^ik 


= 


s G^(x), 


Gfc(x) > 


^ik 


= 


< G.(x), 


Gj^(x) > 


Nijk 


= 


<G.(x)G^(x), Gj^ 


^k 


= 


<f , Gj^(: 


.)> 


Then, 


Eqn. 22 


may be rewritten < 


tv 


+ s Cii^ 


+ MiiJ 


V. 

1 


m 


m 






z: 


z: 


"2 Nijk ''i 


^ - Qi, 
J k 


i=l 


j=i 


J 


k = 1,2 


,3, .... 


m 





(23) 



(24) 



This is a system of nonlinear algebraic equations, where the 
unknown displacements must also satisfy the imposed boundary 
conditions. 
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B. ELEMENT SHAPE FUNCTIONS 



Consider a beam element. There are two nodal coordinates, the 
deflection and slope at each end, and therefore each element must have 
a- total of four degrees of freedom. On the element level these coor- 
dinates are numbered as follows: 

Coordinate 1 = V(o) Coordinate 3 = V(i) 

I I 

Coordinate 2 = V (o) Coordinate 4 = V (t) 

The compatible element shape functions must be cubics and it is 

evident that four independent shape functions are required [Refs. 14, 16]. 
e 

Let G. (x), the element shape functions, be defined as follows, where the 
superscript "e'* refers to the beam element: 




G® = 1 - 3x^ + 2x^ 

l2 ^ 

=x^ -2x^+x (25) 

iZ 1 - 




Note that the function G.(x) represents the variation of V along the 
element in such a way that G|(o) = G^^o) = ^^(i) = G®(1) = 1 and that 
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e e e 0 » 

other evaluations of G (o), G. (1), G. (o), G. (1) vanish; cf. the sketches 

i 1 1 1 

in figure 2 , where x is the local coordinate from left to right and 1 
denotes the length of the beam element. 



C. FORMATION OF ELEMENT MATRICES 



Once the element shape functions are defined, the components of 

the element matrices K^. , Qf, Mf, , Q? and N^ can be computed 

ik ik ik k ^ 

through ordinary integration of polynomials. The following results, 
associated with the linear terms, have been given in many text books 
[Refs. 21 , 14, 25]. 



K' 



ik 





1 * 2 . 6 . - 

l3 j 2 


12 . 

l3 


6 . 


= EI^ 


4. 

-1 


- 6 . 

l 2 


2 . 

-1 




Symmetric 


12 . 

l3 


- 6 . 








4. 

1 ' 


average El) over each 


element. 




1.2 0.1 
1 


- 1 . 2 
1 


C 



e y V 

C =<sG.(x),Gv (x)^= 

ik > -L 't ^ 



. 1333331 -0. 1 
symmetric 



1 . 2 

I 



0.1 

-0. 0333331 

_ 0.1 

0.133333 1 



where s is a constant (the average) over each element. 
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< G®(x), G®(x)> = 



2 2 
. 3714291 : .0523811 . 1285711 : -. 030952 - 1 - 

3 2 

.0095241 .0309521 -. 0071431 ^ 



symmetric 



.3714291 -.0523811 



. 0095241 






.51 

2 

. 0833331 
. 51 

2 

-.0833331 



The components of associated with the nonlinear response are 

not in the literature. They are as follows: 



N.. = 

ijk 


< 


(g (x) G.(x 
^ ^ J 


Kn 

e 


= 


. 3071431 

2 


^211 


- 


. 0384921 


311 


= 


.0642861 


^^411 


= 


2 

-.0170601 


N® 

221 


= 


. 0063491 ^ 


^321 




. 0138891 ^ 




= 


-. 0035711 ^ 


^331 


= 


. 0642861 


N® 

431 


= 


2 

-.0138891 
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_ 


. 0031751^ 


441 


n!-_ 


- 


. 0011911 




^322 


= 


. 0031751 


N® 

422 


= 


-.0007941 




= 


2 

. 0170641 


N432 


- 


-.0035721 


e 




4 


N442 




. 0007941 


n' 


- 


. 3071431 


333 


^433 


= 


-.0384921' 


N® 

443 


= 


. 0063491^ 


e 


N444 




-.0011901 



Notice that, according to Eqn. 23 , the N. coefficients are equal 

I J 

e 0 0 

for any permutation of i, j, k (i. e. , Nu2 = ^121 ~ ^211> ) 

D. FORMATION OF SYSTEM MATRICES 



The system matrices are obtained from the element matrices as 
follows ; 

1 . A correspondence table between local and global coordinates 
is formed 

2 . The element coefficients of any matrix, say the element 

e 

stiffness coefficients are substituted directly into the 

system stiffness matrix Kj.j according to the correspondence 
of element coordinates i, j to system coordinates I, J.(Ref- 14 j 



30 



I 






'¥■ 



t 



i 






0 



V. COMPUTER PROGRAMS 



Equation (13) is originally of the form: 

(EIV")" - (sV)' + kvP -q (29) 

The following block diagram shows how the investigation proceeds 
(Figure 3). Equation 13 or 29 will be solved by the Finite Element 
Method. First consider the case of a classical Winkler foundation, 
for which the theoretical analysis is readily available. The results 
from the Finite Element Method will then be compared to the theoretical 
results. Similarly, the modified Winkler foundation will be solved, also 
by the Finite Element Method. Since there is no theoretical analysis 
for the case of these nonlinear problems, some of the results from the 
Finite Element Method will be compared to those from the Finite Dif- 
ference Method. Finally, in dealing with continuous foundations, cer- 
tain types of soils will be considered. Results obtained by varying 
flexural rigidity of beams, foundation modulus and applied force will be 
presented in the next chapter. 

A. GENERAL PROGRAM ORGANIZATION 

The program may be divided into two separate categories, depend- 
ing on linearity or nonlinearity. The reason for this distinction is that 
for nonlinear problems, the Brown's iteration technique is applied and 
double precision is required. In contrast, a simple Gaussian elimina- 
tion method may be apx^lied to solve the linear systems. 
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Figure 3. Sequence of the problem investigation 
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Figure 4. General. organization of the program 
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1 



( 3 ) 




Figure 5. Matrix formation (STIFFl is used for classical 
Winkler foundation only, STIFF2 is used for all other cases. ) 
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Figure 6. Actual problem solution: 
(a) for linear problem, 

(b) for nonlinear problem. 
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Figure 7. Finite Difference Method solution 
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Each type of program performs four major functions: 

1. Storing all information needed for subsequent problems 

2. Supplying necessary data to each particular problem 

3. a. Forming the element matrices 
b. Forming the system matrices 

4. Solving the algebraic problem. 

This organization is shown in Figures 4 through 6. 

The Finite Difference Method is applied to some problems to check 
the results. Its general flow chart is shown in Figure 7. 

B. SUBROUTINES 

The MAIN program is a control segment that serves to call other 
subroutines in proper order as required. No calculations are performed 
in the MAIN program (See Appendix D). 

1. Subroutine STORE 

This subroutine reads all data available for the problem under 
investigation, including the maximum number of elements, the number 
of nodes, the number of boundary conditions, the number of degrees of 
freedom, the correspondence table of global and local nodes, load 
conditions, beam properties and foundation characteristics. 

2. Subroutine TEST 

This subroutine is used particularly for testing the convergence 
of the solution. In doing so, it selects special information and supplies 
it to the MAIN program which, keeping all information constant, doubles 
the number of elements from one run to another. 
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3. Subroutine CURFIT 



For the case of a nonlinear elastic foundation for which the 
reaction r(x) = kV^ , where p is not an integer, a direct application 
o~f the finite element method via the Galerkin procedure is not efficient. 
Using a least square best fit, this subroutine replaces kV^ by a 
second order polynomial. 

4. Subroutine INCHK 

Depending on various conditions in changing beam flexural 
rigidity and foundation modulus, this subroutine selects information 
from subroutine STORE and supplies it to the MAIN program. This 
subroutine also prints the appropriate echo check. 

5 . Subroutines STIFFl and STIFF2 

Subroutine STIFFl is used for the linear problem and uses 
single precision; subroutine STIFF2 is used for the nonlinear problem 
and uses double precision. Both of these subroutines form the element 
stiffness matrices and assembles the system stiffness matrix. 

6. Subroutine LOAD 

This subroutine forms the element load vector and assembles 
the system load vector. It can be used for both linear or nonlinear 
problems, provided it is in double precision for the later case. 

7. Subroutines BOUNDl and BOUND2 

These subroutines apply the boundary conditions to the system 
force vector and the system stiffness matrix produced by subroutines 
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LOAD, STIFFl or STIFF2, and transform them into an appropriate 
form ready to be solved. Subroutine BOUNDl is used for linear problems 
and subroutine BOUND2 for nonlinear problems. 

8. Subroutine SOLVE 

This subroutine employs the Gaussian elimination method to 
solve the linear system of algebraic equations for the linear problem. 

9. Subroutine RESULT 

This subroutine prints any output data from either linear or 
nonlinear problems. For nonlinear problems, double precision must 
be specified. 

10. Subroutine DIFSOL 

This subroutine solves the nonlinear problem by the Finite 
Difference Method of which the development will be presented in 
Chapter V, part D. This solution is optional and was employed in a 
few cases in order to check the results of the nonlinear problem given 
by the Finite Element Method. 

Subroutine FSOIL 

Similar to subroutine CURFIT, using a least square fit, this 
subroutine replaces the total settlement curve of a natural soil given 
by experimental data [Ref. 30], by a polynomial of second degree. The 
total settlement curve of the sea floor sediment is supposed to have the 
same shape as that of the inland soil but its magnitude is reduced by a 
certain scale. 
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12. Subroutine VGUESS 



This subroutine is designed to give an initial estimate vector 
to the nonlinear algebraic system before any iteration process is made. 
It is desirable that if the first estimate has failed to make the iteration 
convergent, another estimate is provided. In this subroutine, the 
initial estimate must be taken in the form of a second order polynomial. 

13. Function AUX 

This is a nonlinear algebraic equation system translated from 
equation (24), which has to be solved and which serves as an external 
input to subroutine ZSYSTM. 

14. Functions DIFF4 and DIFF8 

These are nonlinear algebraic equation systems translated 
from Eqns. 27 and 28, respectively, which have to be solved by the 
Finite Difference Method and which serve as external inputs to sub- 
routine DIFSOL. 

15. Subroutine LSQPL2 

This library subroutine, from the Naval Postgraduate School, 
Monterey, California, employs a least square best fit method to replace 
any function by a polynomial. 

16. Subroutine ZSYSTM 

This library subroutine, from IMSL, solves nonlinear 
algebraic equations employing the Brown’s method of iteration, which 
is a modified Newton's method and is developed in Reference 13. 
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C. TESTS FOR CONVERGENCE 



The test problems employed here are designed to verify the opera- 
tion of the computer program and test the convergence of the equation 
system. The two categories of test performed were linear and non- 
linear problems. 

An analytic solution and the Finite Difference Method solution were 
used to verify the Finite Element Method solution for linear and non- 
linear cases respectively. The following results (Tables 1, 2, 3) show 
that for an uniform beam of 100 inches in length, with a uniform load, 
the 8-element model gives sufficiently accurate results. These analyses 
show that as the number of elements is increased the finite element 
method results approach the theoretical results. The small difference 
between eight and sixteen element results (less than 1%) suggests the 
solution has converged sufficiently for engineering purposes. For 
symmetric problems (i. e. , symmetric El, k, s and q), advantage 
may be taken of symmetry. In this case, only half the system need be 
considered, thus greatly reducing computer effort. This reduction 
is especially beneficial for nonlinear problems. The time consumption 
may increase from 10 to 30 fold in going from a 4-element model to an 
8-element model, depending on how good the initial guess is. A brief 
description of the test program is shown in Figure 8. 
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Figure 8. Flow chart of a test program 
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0. loJ. 

Deflection 
in inches 

A Winkler 

□ modified Winkler 

^ Winkler 
O modified Winkler 



[ El = S.xlO^ Ibs-in^ 
; q = 22. 5 lbs /in 

) El = 3. xio”^ Ibs-in^ 
J q = 22. 5 Ibs/in 



Deflection 
in inches 



Figure 9. Compared deflection curves of a 
typical inland soil from Figure 15. 
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Deflection 
in inches 




Figure 10. Deflection of a beam on a sea floor model foundation 
compared to that of Winkler foundation. 

(Maximum difference is about 6%) 
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Table 1. Convergence of the Linear System 
(EIV)" + kV = q 

Uniform Simply Supported Beam with Uniform Load 
L = 100 in, El = 3 X 10^ Ib-in'^, 
q = 1 Ib/in, k = 1 Ib/in^, 

Deflection Given in Inches ( /3L = 0. 96) 



Location 


Finite 


Element Method Solution 


Theoretical 


x/L 


Number of Elements Number of Elements 


Solution 




8 


16 




0 . 0 


0.0 


0 . 0 


0 . 0 


0. 125 


0. 016296 


0. 016310 


0. 016301 


0. 250 


0. 029897 


0. 029919 


0. 029906 


0. 375 


0. 038839 


0. 038863 


0. 038880 


0. 500 


0. 41950 


0. 041975 


0. 042178 



Table 2. Convergence of the Nonlinear System 
(EIV")" + kV^ = q 

Uniform Sim ply Supported Beam with Uniform Load 
L =100 in, El = 3 X lo’^ rblin^ 

q = 0. 01 Ib/in, k = 1 Ib/in^ 

Deflection Given in Inches ( /3L = 0. 96) 

Finite Element Method Solution: 



Location 


Number of Elements 


Number of Elements 


x/L 


4 


8 


0. 0 


0. 0 


0. 0 


0. 250 


0. 1451D-03 


0. 0456D-03 


0. 500 


0. 1483D-03 


0. 1485D-03 
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1 1 2 

Table 3. Solution of (EIV>') -f- kV =q 

F. E. M. Compared to F. D. M. Uniform 
Simply Supported Beam with Uniforrn Load 
L = 100 in, El = 6. 123x10^ Ib-in^, 
q = 2.25 Ib/in, Number of Element = 8 



Location 


k = 100 


^L = . 80 


x/L 


F. E. M. 


F. D. M. 




ITMAX = 3 


ITMAX = 3 


o 

• 

o 


0. 0 


0. 0 


0. 125 


0. 185781D-03 


0. 188397D-03 


0. 250 


0. 340909D-03 


0. 345395D-03 


0. 375 


0. 442958D-03 


0. 448565D-03 


0. 500 


0. 478469D-03 


0. 484450D-03 



Location 


k = 500 


^L = 1. 20 


x/L 


F. E. M. 


F. D. M. 




ITMAX = 8 


ITMAX = 3 


0. 0 


0.0 


0. 0 


0. 125 


0. 185776D-03 


0. 188392D-03 


0. 250 


0. 340900D-03 


0. 345385D-03 


0. 375 


0. 442946D-03 


0. 448552D-03 


0. 500 


0. 478456D-03 


0. 484436D-03 



Location 


k = 1000 


^L = 1. 42 


x/L 


F. E. M. 


F. D. M. 




ITMAX = 4 


ITMAX = 3 


0. 0 


0.0 


'0.0' 


0. 125 


0. 185769D-03 


0. 188385D-03 


0. 250 


0. 340889D-03 


0. 345373D-03 


0. 375 


0. 442931D-03 


0. 448536D-03 


0. 500 


0. 478440D-03 


0. 484419D-03 
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Location 

x/L 



k = 2000 

F. E.M. 

ITMAX = 4 



^L = 1. 69 
F. D. M. 
ITMAX = 3 



0. 0 


0. 0 


0. 0 


0. 125 


0. 185757D-03 


0. 188372D.03 


0. 250 


0. 340865D-03 


0. 345349D-03 


0. 375 


0.442900D-03 


0. 448504D-03 


0. 500 


0. 478407D-03 


0. 484384D-03 



Location 


k = 3000 


^L = 0. 87 


x/L 


F. E.M. 


F. D. M. 




ITMAX = 4 


ITMAX = 4 


0. 0 


0. 0 


0. 0 


0. 125 


0. 185745D-03 


0. 188359D-03 


0. 250 


0. 340842D-03 


0. 345324D-03 


0. 375 


0. 442870D-03 


0. 448472D-03 


0. 500 


0. 478374D.03 


0. 484350D-03 
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D. THE FINITE DIFFERENCE METHOD 



Since there is no analytic solution to the nonlinear case, the Finite 
Difference Method was used in this thesis to verify the results of the 
Finite Element Method. The development is restricted to uniform El 
and k. The finite difference approximations are developed in Ref. 15. 
The fourth order derivative has the central difference computational 
molecule as follows: 

(iJL). = J_ [V(j-2) - 4V(j-l)+6V(j).4V(j + l)+V(j+2)] +0(1^) 
where 1 is an element length 

Consider the differential equation of an uniform simply supported 
beam with uniform load on a uniform nonlinear foundation: 

II 2 

EIV + kV = q 

1 1 2 
or EIv" + (EIV) = q 

« k 

Let V = EIV and > then 

the above equation may be rewritten as follows: 

V" + A V ^ - q = 0 



A computational molecule model applied at node j yields: 

_i_ v"(j-2) _ 4V^^j-l) + 6V*(j) . 4V*(j+l)+V*(j+2) 

X4 

+ A V*2(j) _ q(j) = 0 - - 



or: 



V (j-2)-4V (j.l) + 6V (j) : 4V (j+1) + V (j+2) 
4 ^'2 4 

+ A X^V (j) - X q(j) = 0 



(26) 
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where j ranges from 1 to 2 for 4-element model and from 1 to 4 for 
8-element model. 



A. FOUR- ELEMENT MODEL 



V(-1) = -V(l) V(0) 



V(l) 



V(2) 



V(l) 



• — — 



X 



6 V(l) - 4 V(2) + A X^V^(l) _ qX'^ = 0 
-8 V(l) + 6 V(2) + A X'^V^'(2) _ qX"* = 0 



B. EIGHT- ELEMENT MODEL 



V(0) 

# 



(27) 



-V(l) V(0) V(l) V(2) V(3) V(4) V(3) V(2) V(l) V(0) 



(28) 



" X ' 

5V(1) - 4 V(2) + V(3) + A x"^V^(l) - qx"^ = 0 
-4V(1) + 6V(2) - 4V(3) + V(4) + A X^V^(2) - qX^ = 0 
V(l) - 4V(2) + 7V(3) - 4V(4) + A X^V^(3) - qX^ = 0 
2V(2) - 8V(3) + 6V{4) + A X^V^(4) - qx"^ = 0 



Notice that the superscript has been omitted for simplicity; the 
result V(j) must be multiplied by (EI)”^ to get the actual deflection. 

The above development is restricted to a uniform beam on a founda- 
tion with uniform k. In this case the finite difference method has 
distinct advantages over the finite element method. In the case of 
variable El, k, and s the finite difference method is not easily developed 
and the advantage shifts to the finite element method. 
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VI. ILLUSTRATIVE PROBLEMS AND RESULTS 



This chapter presents the results of two investigations. 

The first investigation attempts to compare the behavior of three 
foundation models, the Winkler, the modified Winkler and the mixed 
mode foundations. The comparison is limited to the case of a simply 
supported uniform beam on foundations with uniform k and s properties. 
Table 4 gives values for the parameters considered in this study. The 
value of ^L varied from 1. 62 to 3. 98 i. e. , beams of medium length. 
The results obtained apply only to the particular simply supported 
systems considered. Other problems may yield other results. 

The second investigation seeks to establish that the Finite Element 
Program developed in this thesis can solve problems with variable El, 
k, s and q. As a sample problem, the case of a free-free beam with 
variable stiffness, variable foundation modulus, and variable load was 
considered. The jSL for this illustrative problem was 3. 5. 

A, UNIFORM SIMPLY SUPPORTED BEAMS WITH UNIFORM LOADS 
ON SOIL FOUNDATIONS 

Figures 9 and 10 present the deflection of a simply supported beam 
on soil foundations which are either inland soil or sea floor sediment. 
The range of parameters considered is shown in Table 4 below. 
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Table 4. Information Data 



p 


Loading 

Conditions 

q 

(Ib/in) 


Beam 

Properties 

El 

(Ib-in^) 


Foundation Characteristics 


Winkler 

type 

k( 


Continuous 

type 


Mixed mode 
(natural soil) 


2 


1 


3x10^ 


1 


k 


s 


inland 


1. 6 


2.25 


3x10® 

Q 


50 


Ib/in^ 


lb 


soil and 


1.2 


22. 5 


6x10 


100 






sea floor 


1. 


225. 


9x10® 


500 


8. 6 


18. 3 


were 


. 9 






1000 


500 


1000 


chosen 


. 85 






1500 






as 


. 8 






2000 


3000 


6000 


examples 








2500 














3000 
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Appendix B contained the tabulated results of Figures 9 and 10 as 
well as data on the two foundations, inland soil and sea sediment. The 
results show that the difference between the Winkler hypothesis and 
other hypotheses is negligible for simply supported beams with uniform 
properties and loads. This may not be true for beams with other 
boundary conditions. 

B. A FREE-FREE NON-UNIFORM BEAM PROBLEM 

To show how the Finite Element Program can be used for a general 
case we consider the following problem: 



lb 




2 



5000. 




4000. Ib/i 



2000 
lb /in 




400 

Ib/in^ 




k 



2 



i 



50. in 
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A symmetric problem was considered only to minimize the 
computer effort. The result is given by Table 5 and Figure 11. 
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Table 5. Deflection of a 
Free-Free Non-Uniform Beam 



Location # 1 




Location # 
1 
2 

3 

4 

5 



Deflection, inches 
1. 19150 
2, 09179 
2. 87908 
3. 40052 
3. 58134 
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Deflection, 

inches 



0 . 



1 . 



2 . 



3. 




4. 



1 



2 3 



4 



5 Location # 



Figure 11. Deflection of a free-free 
non-uniform beam 
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VII. CONCLUSIONS AND RECOMMENDATIONS 



The computer program presented in this thesis provides an accurate 
and reliable means for solving a variety of one-dimensional problems of 
beams on nonlinear foundations. The use of this program and efforts to 
increase its versatility are encouraged. 



A. REMARKS ON THE PROGRAM 

The following observations have been made during the development 
of the program and the investigation of various problems 
1 . Curve Fitting 

It was assumed that the reaction of the foundation is proportional 

to , where V is the deflection and p any positive real number. 

Since the Finite Element Method via the Galerkin procedure is not easily 

applied directly for the case when p is not an interger, has been 

2 

replaced by a least square fit polynomial of the form b^ + b2V + ^3^ • 

The reaction of the foundation must be zero when there is no deflection, 
therefore the greatest error may come from the non-zero residual 
constant b^ (Figure 12). Suppose the least square fit curve 
b^V + ^3^ intercepts the curve r^ = at the point V = V”. If the 

deflection of the foundation falls in the vicinity of V, the result is highly 
accurate. The result is less accurate, if the deflection is much larger 
or smaller than V . One must estimate how large or small the deflection 
will be, in order to select a ’’best fit” curve. 
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Figure 12. Accuracy of curve fitting 

2. Initial Estimate of Deflection 

Iterative procedures are among the more popular schemes for 
solving systems of nonlinear algebraic equations. Such methods require 
an initial estimate of the solution. In this work the initial estimate was 
taken in the form of a second order polynomial in the subroutine 
VGUESS. Not all initial estimates will yield a solution; therefore it is 
desirable that if one initial estimate fails, it be replaced by another 
one. This is accomplished in subroutine VGUESS which changes the 
scale of the initial estimate function. 

3. Computational Considerations 

In any computer effort there are two major concerns, the 
space required for a program and the time necessary for solution. 
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For the linear Finite Element problem the storage area required is 

proportional to m X M, where m is the number of coordinates and M 

is the bandwidth of the system. A change in m or M yields a correspond. 

ing linear change in the size of the computer storage. In addition the 

time for solution (for a Gauss elimination type scheme) is of the order 

2 

of m X M for a symmetric system. These considerations show the 

disadvantages which result when m or M increase, and are well known. 

In the case of nonlinear problems the effects of increasing m 

and/or M on computational effort are even more pronounced. For 

example the space requirement associated with the nonlinear terms is 

2 

m X M. Hence an increase in m from m to m increases the storage 

12 ^ 

m Z 

in the ratio ( . This is times greater than in the linear case. 

mj mj 

Moreover, the computational effort is greatly affected since a large 
number of terms greatly increases the number of iterations. For 
example in the solution of an 8- element model (i. e. , m = 18, M = 6) 
non linear problem, the storage for the array is about 8000 bits (single 
precision) and the computer time was about 60 seconds. When the same 
problem was modeled by a 4-.element model (i. e. , m = 10, M = 6), the 
storage was about Z400 bits (single precision) and the computer time 
was about 6 seconds; a reduction in storage in the ratio of 3 to 1 and a 
reduction in computational time in the ratio of 10 to 1. 
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B. REMARKS ON THE WINKLER HYPOTHESIS 



Let r = kV^ be the total settlement curve, i. e. the deflection 
n 

versus reaction curve of a real foundation, which can be obtained by an 

experiment [Ref, 30]. Let the straight line = kV, where k is the 

foundation modulus, represent the behavior of the Winkler foundation 

(Figure 13). Let V correspond to the point at which the two curves r^ 

and r intersect. Consider the area under the curve r = kV^ 
n ^ 

between 0 and Vg, This area is equal to: 



rv. 



Nr 



rdV 



kV^ dV = 



P + 1 



V 



P+1 



(30) 



Consider the functional: 



■D 



•2 



(EIV) + (° V) + 



P+1 



V 



P+1 



- qV)Jdx 



(31) 



where L is the length of the beam. The general field equation (29) is 

the Euler *s equation of the above functional. Therefore, to solve the 

equation (29) is the same as to minimize the functional (31), This 

functional represents the total potential energy of the system, including 

the beam, the foundation and the applied load. 

The third term under the integral is the energy due to the reaction 

of the foundation. According to Eqn, 30, the energy (31) is a function 

of the area under either the r. or r curve, depending on which founda- 

1 n 

tion is being used. Neglecting the effect due to the shear force in the 
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r 




Winkler foundation 



r 

lbs /in 



nonlinear foundation 




V, inches 



Figure 13. Comparison of a Winkler foundation 
to other types of foundations 



foundation for a while, for a specified beam and a given applied load, 
one can see that if the deflection is small, say less than V, the Winkler 
foundation deflects more than the nonlinear foundation does for equal 
supporting effort. Now, if the effect of the shear force in the foundation, 
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i. e. , — V , is taken into account, one can expect that the deflection 

12 

should be less. The shear force of the foundation depends on (sV) 

»‘2 

which is usually much smaller than (EIV) and may be neglected. This 
discussion agrees with the results in two examples of two different 
types of soils given in Chapter VI, part A. 



C. RECOMMENDED FUTURE STUDIES 

1. Following a similar procedure, the problem of plates or shells 
on two- or three- dimensional nonlinear foundations may be investigated. 

2. This thesis dealt with nonlinear foundations which are either 
continuous or discontinuous. Real foundations may be some combination 
of these foundations; hence there is a need to bridge the gap and consider 
foundations where the deformation is partly localized and partly 
continuous. 

3. In the present investigation the nonlinear term kV^ was 

2 

approximated by the second order polynomial, bj + b 2 V + . This 

restriction results from the finite element approximation of V(x) by 
m 



E 



V.G.(x). 
1 1 



Future studies might consider ways to treat kV^ or in 



fact the more general case of an arbitrary nonlinear function of V, 
directly. 



61 



p 

- 

I 




APPENDIX A 



CONTINUOUS FOUNDATION ANALYSIS 



A. GENERAL FORMULATION 



This formulation is limited to plane stress theory only. Let u(x, y) 

and V (x, y) be the displacements at a certain point of the foundation, in 

the X- and y- directions.’ In general, an approximate solution can be 

obtained by expanding u(x, y) and v(x, y) in a finite series: 
m 

u(x, y) ^ u^(x) (p^(y) , 

i=l 



v(x, y) 



= 23 V (x) 

j = l 



ipAy) 



(32) 



where the dimensionless functions ^.(y) stnd ^j(y) are known and 
the functions U^(x) and Vj(x), which have the same dimensions as 
u(x,y) and v(x,y), are unknown. 

From the theory of elasticity, in the two-dimensional plane stress 
case, the stresses and strains are related as follows [Ref. 3]: 



'XX 






^ ^xx ^ ‘'f^yy - ^ 



-yy 



f 

1 -^ 



( H 



+ ) 



(33) 



f ^xx ^yy 



xy 



yx 



2(1 + 



xy 
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XX 



•yy 






xy = 



du 

c)x 

^y 

*^yx = 



^y dx 



\ 



> 



( 34 ) 



where and are the modulus of elasticity and the Poisson*s ratio of 

the foundation, e and ^ are the strains in the x- and y- direction 
XX yy 

and is the shearing strain. 

In terms of the series (32), the relations (33) may be rewritten: 




ayy = 



Ef 



A - 1 



I 



u . <p 







(35) 



’’xy = Tyx 




2(1+ l^£) 




] 



where the ''prime** denotes the derivative of the function with respect 
to its own argument. 

4 

From the foundation model shown in Figure 1, consider now an 
elementary strip of length Ax (Figure 14). The functions U^(x) and 
V.(x) can be obtained from the equilibrium conditions. According to 
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z 



Ax 




Figure 14. Force system on a strip 
of the foundation 
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4 




I 







Lagrange's principle of virtual displacements, ^ the total work of all 
internal and external forces acting on this strip due to any virtual dis- 
placement is zero. Let: 
u (x, y) = 



V (x, y) 






(36) 

(37) 



be the virtual displacements, where and Vj are arbitrary constants 
at that location x of the strip. 

Since there is no force applied on the surface of the foundation, 
and the bottom of the foundation rests on a rigid base, there is no work 
done on these two faces. It is assumed that all properties of the founda- 
tion remain constant in the z-direction, hence there is no pressure 
gradient in the z-direction and therefore the total work done on two 
faces perpendicular to the z-direction is zero. The system of forces 
which must be taken into account is shown in Figure 14. 

The external forces are the result of normal stresses 






9a 



°xx + ax 



^ Ax, the shearing stresses T"xy> ^xy 



Q'Txy Ax and the distributed forces, X (x, y) in the x-direction and 
9x 

Y(x,y) in the Y-direction. 



The internal forces are the result- of normal stresses 



TY 



and 



shearing stresses work done by these internal forces is 



A virtual displacement is any displacement which is consistent 
with the forces and constraints [Ref. 2]. 
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the strain energy of the strip [Refs. 4, 9]. 



The total work done on the strip by any m + n virtual displacements 
is then given by the following expression; 



I 

/ 

0 

r 

•'n 



^x 



A X ) (wdy) - 



+ I X (x, y) ( A xwdy) + 






I A X ) (wdy) Vj 



(38) 



H 



•'O X 

H 



^ ^yy ^ yy^ ^ ^ 



0 



I 



Y (x, y) ( Ax wdy ) Vj = 0 



(x l|^y3^...y m , J l|2^3j...^n) 
where 



y 



*<)Ui 






xy 



U. <P (y) 



(39) 



t) V 



:yy 






f. 



H 






XX 



■2>y 

Ax) (wdy) u. 



V. 0 (y) 



(40) 



= work of external forces ( ^ ^^Ay )(wdy); 

‘«>x 



xy 



I ( ^ xy Y^y) { Ax wdy) = strain energy due to shearing stress ^xy» 
^0 - - 

fi2l: 

■Jq X 

_ 

I ( ^ ) { Ax wdy) = strain energy due to normal stress 

■^0 ^yy yy yy 



Ax) (wdy) = work of external forces (*^^xy Ax ) (wdy) 

^ *b 



X 



fx., 

J 0 



yX Ax wdy) 



u.. r“ 

‘ J ' 
-'0 



Y(x, y) ( Ax wdy) Vj ~ work of body forces. 
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After dividing through oy wAx and substituting tr^, , y^y> ~^yY 
(36), (37), (39) and (40), the expression (38) may be rewritten as follows: 



m 



i=l 









1 ®xx _ 

1 ^i(y) dy - 

-'o ^ X 

4 - 


Txy i (y) dy + . 


1 

) 0 X (x, y) U i(y)dy 




_ 



+ ^ \ f v; (y)dy - J H^^Y) ^y) 

j=lL^ JJ y Q 



^ Y (x, y)Vj v|/ j (y) d>j= 0 



(41) 



or 



m 

Z 

i=l 



■u . [ 

1 



rii 1 

1 1> <^xx (y) dy - \ “^xyCp^ (y) dy + \ X (x, y )^ (y )dy]> 

-^0 T)x ‘ J 0 *'o ^ 



+ Z I [ I„ lljs >)- j 

j=l ' J ° -Dx 



(y) dY - 



f ^ ' r ) 

Jq <Tyy 'f j “'o Y(x, Y)tj(Y)dY]| 



(42) 



Since U and Vj are arbitrary, the above expression leads to the two 
equations in the x- and y-directions : 



\ " ^(Txx cp^(y)dy - \ Txyy.(y)dy+ C X (x, y ) (^. (y)dy = 0 (43) 

*^0 -^0 ^ Jo 



r ?>rxy 

•'n 



-H 



rn , 

j CTyy v^^(y)dy + J ^ Y(x, y) j(y)dy = 0 (44) 



0 






;(Y)dy 
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I 




Substitution of (35) into (41) and (42) leads to a system of m + n 



second order differential equations in U^(x) and follows: 



^ E 

c) 


m 

i=i 


m 

j=l 




m 


n 


.{ ‘ 




+i:^ 


Jq 2(1+'J'£) 


i=l 


j=i 






H 



1 ^. dy = 0 



>(45) 



•H 






n I 



[ V f i +51' f j 

/(.-Sx 2(lWf) -^1 



1 = 1 j=l 



H 



dy + Yv|/^ dy 



= 0 



y 



(k = 1, 2, 3, . • . m 
t = 1,2,3,.. .,n) 

These are the two governing equations of a two-dimensional 
foundation which has a depth H. 



B. SPECIAL CASE 

For the limited scope of this thesis, let the horizontal displace- 
ments u be negligible. This means that any virtual displacements in 
the x-direction are also negligible. Since TT are arbitrary constants, 
Eqn. 36 implies that ^^(y) must be identically zero. Only the last 
equation of the system (45) remains and it reduces to: 
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^ . E f TT 14 » 

i=l 



2(1 + Vf) Jo 

H 



Tj 



.- dy 



1 



+ j Y (x, y) v|/ j dy = 0 



(46) 



where j = 1, 2, 3, . . . , n, and E£ and are assumed to be constant. 
The above expression may be rewritten: 



m 



E 



i = l 



•a 1( 






or 



7)x 2(1+Vf) 

r 



r = 
J 



\ yitj Vj ] . [ Ef j dy] V. 

0 1 - V ^ *^0 



Y (x, y) j dy = 0 



(47) 



m 



E 



• f 



i=i 



'(^ijV;) . kyV; 



(4 8) 



where 



•H 



= 



r 

I Y (x,y) dy 

'^0 



(49) 



is the reaction of the foundation. 



«ij = 



= 



^f 


1 H'i Tj dy 


2(1 + Vf) 


Ef 


1 * f 


1- vj 


J S' i f j 
•'0 



(50) 



( 51 ) 



69 



The coefficients Ibs/in, and k— , Ibs/in , are the two characteristics 

of a foundation [Ref, 5], The former accounts for the shearing stress 
distribution in the foundation material, and the latter is equivalent to the 
Wihkler constant of proportionality. They are however, dependent on 
how the 0. functions are chosen and are consequently not independent 
of each other. 

For simplicity, assume that the foundation has an infinite depth, 
which is usually true for any soil layer compared to the deflection of 



most structural systems. Let the one-dimensional function ^(y) 
satisfy the boundary conditions: 

^(o) = 1. (52) 

^{oo) = 0. (53) 

Then the displacement at any point of the foundation is: 

v(x, y) = V (x) ij) (y) (54) 



in which V(x) is the deflection at the surface of the foundation, where 
y = 0. Furthermore, assume that the displacements decrease ex- 
ponentially with the depth of the foundation, then may be 

expressed as: 

- Xy 

ip(y) = e (55) 

where X is a known constant which depends on the properties of the 
foundation. The values of s and k* are then given by (50) and (51): 



E, 



S = 



^ CO 



2(1 + Vf) 



2X 



^dy = 



E, 



4 X( 1 + V£) 



Ibs/in (56) 
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k*= (■“’ (e'^'')'^dy= lbs/. 3 (57) 

Equation (48) may be rewritten 

r = (sV’)' - (58) 

where s and l^are given by (56) and (57). 



C. NUMERICAL APPLICATIONS 

In accordance with soil mechanics, taking the length of a loaded 
area as unity, Reference 7 gives an empirical relation between bearing 
load and settlement of a foundation as follows: 

1 E, 

r = 1 L_ V (59) 

M 1- Vf 

where pi is called the influence coefficient, which depends on the 
shape of loaded area and the properties of the foundation. The values 
of pi are given in Table 6. [Ref. 6]. Comparison of the Winkler 
hypothesis 

r = kV (60) 

to equation (59) yields: 

k = _J_ ^f (61) 

M 1-V^f 

Eqns. 57 and 61 show the similarity of the theoretical analysis and 
the empirical result. The value X can be computed by equating 
(57) to (61): 

X = 2pi (62) 
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Table 6 

Influence Coefficients 



{After Z. Wilun and K. Starzewski [Ref. 7]) 



Shape of 
foundation 


Flexible foundation 


Rigid foundation 


Settlement 
of centre 


Settlement 
of corner 


Mean 

Settlement 


Settlement 


Mo 


Me 


M 


Mr 


Circular 


1. 00 


0. 64 point on 
perimeter 


0. 85 


0.79 


Square 


1. 12 




0.95 


0. 88 


Rectangle 










2 X 1 


1. 53 




1.30 


1. 22 


5 X 1 


2. 10 




1. 83 


1.72 


10 X 1 


2.53 




2.25 


2. 12 


100 X 1 


4.00 




3.69 


-- 
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For the purpose of this thesis, consider a loaded area of shape 
wx A X, where Ax is taken as a unit length, one inch, and w is 
equal to the width of the foundation (Figure 14). Assume that w 
ranges between 5 inches and 10 inches. From Table 6 , 1. 83 < ^ < 2.25. 
Let p = 2 , then, from (62) \ = 4 , In order to determine s and k, 

the modulus of elasticity and the Poisson*s ratio of the foundation must 
be known. 

The modulus of elasticity of rocks ranges from 100 to 10,000 
2 2 2 

mn/m (7mn/m = Iklbf/in ) and the Poisson*s ratio from 0. 1 to 0. 3 
[Refs. 10, 7 ]. Deep sea sediment cores have modulii of elasticity 
ranging from 0 to 12.4 psi [Ref. 11]. The modulii of elasticity of 
beach soils vary from 30 to 570 psi, while those of inland soils vary 
from 3, 000 psi to 11,500 psi [Ref. 7]. Most soils have a Poisson^s 
ratio near to 0. 5 [Ref. 12]. The values of s and k computed from 
(56) and (57) for various types of soils are given in Table 7. 

Table 7 



s and k Values Depending on Various Types 
of Soils (Loaded Area of Shape Ax: - 1", 





5" < 


w < 10") 






Sea Floor 
Sediment 


Beach Soils 


General 
Inland Soils 


Ef, psi 


12. 4 


30-570 


3000-11, 500 




0. 5 


0. 5 


. 45 


s, lbs 


18. 6 


50 


5, 000-21, 000 


k, psi 


8. 3 


20-380 


1800-7200 
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APPENDIX B 



DATA AND RESULTS FOR PROBLEMS OF CHAPTER VI, PART A 

The following example refers to Figure 23-10 of Reference 30. 
Assume that the piles used for this soil have 24-inch diameters. The 
reaction for the foundation is 

r = (Test load) w lbs /in 

irD^74 

With D = 24 inches and w = 6. 5 inches we have the following table. 



Deflection Versus 


Table 8 
Reaction of a 


Natural Soil 


From Experimental Data 




[Ref. 30] 




lbs /in 


Deflection 


Test load 




(inches) 


(tons) 




0. 10 


20. 


574. 6 


0.25 


40. 


1149. 2 


0. 40 


50. 


1436. 5 


0. 60 


60. 


1724. 5 


0. 90 


70. 


2011. 8 


1. 30 


73. 


2097. 6 



A curve of this data is shown in Figure 15. 
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f 




r(kip/ in) 

inland 

soil 

3. 

2 . 

1 .. 

c 

0 , 




Experimental data curve 
D Least square fit 

^ Assumed foundation behavior curve 

Approximate Winkler hypothesis, k — 3000. Ibs/in^ 

Figure 15. Foundation reaction versus deflection 
(After Spangler, M. G. , [Ref. 30]) 
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A. A TYPICAL INLAND SOIL 



Figure 15 shows reactions versus displacements of a typical inland 
soil (read scale on left hand side). 

2 

Assume that the Winkler foundation modulus is about 3000 Ib/in . 

The deflection results, obtained by the Finite Element program, of 

a simply supported uniform beam with uniform load resting on this type 

of soil are given by Table 9. 

Column 1: Solution of (EIV")*' + kV = q 

1 1 

Column 2: Solution of (EIV’*) + r(V) = q where r (V) is replaced 

by a least square fit curve of the reaction versus deflection curve given 
in Figure 15. 

M , I 

Column 3: Solution of (EIV") - (sV ) + r{V) = q where s and k 

corresponding to a particular foundation is given by Table 7. 

B. AN ASSUMED SEA FLOOR MODEL 

Since there is insufficient experimental data of sea sediment, it 

will be assumed that the reaction versus deflection curve of sea 

sediment has the same shape as that of inland soil. According to 

Table 4 and Reference 11, assume that the Winkler foundation modulus 

2 

of sea floor is about 8. 3 Ib/in . It will be assumed that the inland soil 
reaction versus deflection curve is scaled down by 8.3/3000 to obtain 
the sea floor reaction versus deflection curve. (Figure 15. read right 
hand side scale). 
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The deflection results of a simply supported uniform beam with 



uniform load resting on this type of sea foundation model is given by 
Table 10. 
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Table 9. 


Compared Foundation Deflection 




of Various Hypotheses for a Typical Inland Soil 




Represented in Figure 15 




a) 


L = 100 in, El = 


3 X 10^ Ibs-in^, q = 22. 


5 lbs /in 




k = 3000 Ibs/in^ 


, s = 6000 lbs /in 




Location 

x/L 


(1) Winkler 


(2) Modified 
Winkler 


(3) Mixed mode 
foundation 


0. 0 
0. 125 
0.250 
0. 325 
0. 500 


0.0 

0. 003690 
0.006390 
0. 007898 
0. 008370 


0. 0 

0. 350190D-02 
0. 603699D-02 
0. 744817D-02 
0. 788845D-02 


0. 0 

0. 349588D-02 
0. 602663D-02 
0. 743554D-02 
0. 787515D-02 


b) 


L = 100 in, El = 


7 2 

3 xlO Ibs-in , q = 22. 5 


lbs / in 




k = 3000 Ibs/in^ 


, s = 6000 Ibs/in 




Location 

x/L 


(1) Winkler 


(2) Modified 
Winkle r 


(3) Mixed mode 
foundation 


0. 0 
0. 125 
0.250 
0. 375 
0. 500 


0. 0 

0. 005513 
0. 007740 
0. 007988 
0. 008010 


0. 0 

0. 524425D-02 
0. 728582D-02 
0. 749146D-02 
0. 739400D-02 


0. 0 

0. 522121D-02 
0. 726411D-02 
0. 748313D-02 
0. 739234D-02 


c) 


L = 100 in, El = 
k = 3000 Ibs/in^ 


7 2 

: 3 X 10' lbs -in , q = 225 lbs /in 
, s = 6000 lbs /in 


Location 

x/L (1) Winkler 


(2) Modified 
Winkler 


(3) Mixed mode 
foundation 


0. 0 
0. 125 
0. 250 
0. 375 
0. 500 


0. 0 

0. 055125 
0. 077400 
0. 079875 
0. 080100 


0.0 

0. 536074D-01 
0. 749189D-01 
0. 773826D-01 
0. 764925D-01 


0. 0 

0. 533658D-01 
0. 746865D-01 
0. 772876D-01 
0. 764673D-01 
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Table 10. Compared Foundation Deflection 





of Various Hypotheses for a Sea Floor Soil 




Represented by Figure 15 




L = 


7 

100 in, El = 3 X 10 


2 

lbs - in 




' a) q 


= 1 lbs /in 






Location 




Modified 


Mixed mode 


x/L 


Winkler 


Winkler 


foundation 


0.0 


0. 0 


0. 0 


o 

• 

o 


0. 125 


0. 131617E-01 


0. 130091D-01 


0. 130029D-01 


0.250 


0.241054E-01 


0. 238239D-01 


0. 238124D-01 


0. 375 


0. 312727E.01 


0. 309054D-01 


0. 308904D-01 


0. 500 


0. 337603E-01 


0. 333630D.01 


0. 333467D-01 



b) q = 


2.2 5 lbs /in 






Location 

x/L 


Winkler 


Modified 

Winkler 


Mixed mode 
foundation 


0.0 
0. 125 
0.250 
0. 375 
0. 500 


0. 0 

0. 296138E-01 
0. 542371E.01 
0. 703635E-01 
0. 759606E-01 


0. 0 

0.293785D-01 
0. 538043D.01 
0. 698002D-01 
0. 753519D-01 


0. 0 

0. 293644D-01 
0. 537781D-01 
0. 697660D-01 
0. 753150D-01 


c) q = 


4. 5 lbs /in 






Location 

x/L 


Winkler 


Modified 

Winkler 


Mixed mode 
foundation 


0. 0 
0. 125 
0. 250 
0. 375 


0.0 

0. 592276E-01 
0. 108474 
0. 140727 


0. 0 

0. 591543D-01 
0. 108346 
0. 140568 


0. 0 

0. 591253D-01 
0. 108292 
0. 140498 



d) q = 6. 75 Ibs/in 



Location 




Modified 


Mixed mode 


x/L 


Winkler 


Winkler 


foundation 


0 . 0 


0.0 


0. 0 


0. 0 


0. 125 


0. 8884616E-01 


0. 893441D-01 


0. 892998D-01 


0. 250 


0. 162712 


0. 163656 


0. 163574 


0. 375 


0.211090 


0.212344 


0. 212237 


0. 500 


0.227882 


0.229248 


0. 229132 


e) q = 


9 lbs /in 







Location 




Modified 


Mixed mode 


x/L 


Winkler 


Winkle r 


foundation 


0. 0 


0. 0 


0. 0 


0. 0 


0. 125 


0. 118455 


0. 119966 


0. 119906 


0.250 


0. 216949 


0. 219677 


0. 219657 


0. 375 


0. 281454 


0.285174 


0. 285028 


0. 500 


0. 303843 


0. 307884 


0. 307726 



f) q = 


15. 75 lbs /in 






Location 




Modified 


Mixed mode 


x/L 


Winkler 


Winkler 


foundation 


0.0 


0. 0 


0. 0 


0. 0 


0. 125 


0. 207297 


0. 214630 


0. 214517 


0.250 


0. 379660 


0. 393300 


0. 393091 


0. 375 


0. 492545 


0. 510477 


0. 510204 


0. 500 


0. 531724 


0. 551182 


0. 550887 


g) q = 


22. 5 lbs /in 




- 


Location 




Modified 


Mixed mode 


x/L 


Winkler 


Winkler 


foundation 


0. 0 


0.0 


0.0 


0. 0 


0. 125 


0.296138 


0. 313991 


0. 313816 


0.250 


0. 542371 


0. 375550 


0. 575277 


0. 375 


0. 703635 


0. 747221 


0. 746800 


0. 500 


0. 759606 


0. 806886 


0. 806430 



80 



APPENDIX C 



COMPUTER PROGRAM NOMENCLATURE 

A complete listing and description of all variables used in the 
program is not practical. The items listed in this appendix are common 
to several areas of the program and will assist the reader in a study of 
the program. For convenience items are listed by variable type. The 
definition or function and dimension, if applicable, of each item is 
given. 

A. INTEGER CONSTANTS 



KA 


maximum number of data of foundation modulus 


KEI 


maximum number of data of beam flexural rigidity 


KQ 


maximum number of loading conditions 


KS 


maximum number of various type of soils 


KATEST 


specified foundation modulus for test program only 


KETEST 


specified beam flexural rigidity for test program only 


KQTEST 


specified loading condition for test program only 


ITMAX 


maximum number of iterations 


NBC 


number of boundary conditions imposed 


NCODE 


code number, if NCODE is not equal to 1, the shear force 
in the foundation material is not taken into account 
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NDOF 


number of degrees of freedom 


NELEM 


number of elements 


NELMAX 


maximum number of elements 


NNOD 


number of nodal points 


NSIG 


number of significant figures 



B. REAL CONSTANTS 



ALPHA 


foundation modulus k 


El 


beam flexural rigidity 


Q 


applied load 


SSOIL 


shear force in a soil foundation 


TL 


total length of the beam 


X 


beam element length 



C. VECTORS 



ASAV(IO) 


a set of foundation modulii 


B(10) 


coefficients of a least square fit polynomial 


EF(IO) 


element force vector 


EISA V{ 10) 


a set of beam flexural rigidities 


IDBC(20) 


identification number of boundary conditions 


PO(IO) 


a set of power p’s in the term kV 


QSAV(IO) 


a set of applied loads 


SK(IO) 


a set of foundation modulii of natural soils 


SS(IO) 


a set of shear forces of natural soils 
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1 




SYS F(20) 


system force vector 


THETA(20) 


slopes at the nodal points 


V(20) 


deflections at the nodal points 


WA(600) 


working vector 


D. MATRICES 
DC(4, 4) 


- matrix, referred to Eqns, 23 and 24 


DM(4, 4) 


- matrix, referred to Eqns. 23 and 24 


DN(4, 4, 4) 


N— - matrix referred to Eqns. 23 and 24 


EK(4, 4) 


element stiffness matrix referred to Eqns. 23 



and 24 



ICORR(20, 4) 


correspondence table of local and global points 


SYSC( 20, 20,20) 


system matrix, referred to Eqn. 24 


SYSK(20,20) 


system ^ik- combined matrix. 




referred to Eqn. 24. 
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4 






oooooo oooo ooooooo oooooo ooooooooo 



• ,TIME = 1 



//TAM JOB ( 1714, 0733, NS34) , • 
// EXEC fORTCLO, REGI0N=200K 
//FURT.SYSIN DO * 



THIS MAIN PROGRAM MAY BE USED FOR EITHER 
LINEAR OR NONLINEAR PROBLEMS , ACCORDING TO 
APPROPRIATE INSTRUCTION . FOR LINEAR 
PROBLEMS, OMI T ALL 'IMPLICIT REAL=4=8 ( A-H, G-Z ) ' 
CARDS IN THE MAIN AS WELL AS IN ANY SUBROUTINE 



IMPLICIT REALMS! A-H,0-Z) 

COMMON EK(4,4J ,DC(4,4) ,DM(4,4J ,DN(4,4,4) , 

1 IC0RR(20,4) ,SYSK( 20,20) »SYSC( 20,20,20) , 

2 SYSFI 20) ,£F(4) , IDBCI 20) ,B( lOi , 

3 NNOD, NELEM, ND0F,NBC, TL , AAl, AA2, 

4 X , El , Q, ALPHA, SSOI L,A1 , A2 
DIMENSION WAI600) , THETAI20) ,V( 20) , 

1 PO(IO) ,SK(10) tSS(lO) , 

2 ASAVI 10) ,£ISAV( 10) ,QSAV( 10) 

EXTERNAL AUX 

1 FORMAT! 'I* ) 

THE NCODE NUMBER IS USED TO CONTROL THE 
OUTPUT RESULTS DEPENDING ON THE TYPES 
OF PROBLEMS . NCCDE=l FOR LINEAR 
PROBLEMS . NC0DE=2 FOR NONLINEAR PROBLEMS. 

NC0DE=2 

CALL STORE I MDDF, NBC, T L, NELEM, 

1 KATEST,KETEST,KgTEST,KA,KEI, 

2 KQ,KS,EPS,NSIG, 

3 ASAV,QSAV,E ISAV,PO,SS,SK, ICORR) 

THE FOLLOWING PARAMETERS MAY BE USED IN 
THE CASE OF UNIFORM BEAM WITH UNIFORM 
LOAD AND MAY VARY ACCORDINGLY TO THE 
beam FLEXURAL RIGIDITY , LOADING CONDITION 
AND TYPES OF FOUNDATIONS, 

I P=1 
JP = l 
KP = 1 
LF = 1 
PRINT 1 

CALL I NCHK ( NELEM, NNOD, NBC, MDOF, 

1 X , TL, El » 0, I P, JP, KP, I SOIL, 

2 SSOIL, ALPHA, NCGOE,VSCALE, 

3 ASAV,EISAV,QSAV ,SK,SS, IDBC,V, PO) 

FOR THE CASE OF LINEAR PROBLEMS , REPLAC E 
•CALL STIFF2' CARO BY 'CALL STIFFl* CARO 

CALL STIFF2 
CALL LOAD 
CALL BOUND 

FOR THE CASE OF LINEAR PROBLEMS , RE PL ACE 

THE FOLLOWING CARDS BY 

CALL SOLVE(SYSK,SYSF,V,NNOO) 

CALL RESULT (NNOD ,V ,£l , lER, ITMAX,NCOUE) 

NG = 0 

4 CALL VGUESSINNOD, NELEM, VSCALE, ITMAX,V) 

CALL ZSYSTM( AJX,EPS,NSIG,NNOD, V, ITMAX,WA ,PAR,IER) 
IF(IER.EQ.O) GO TO 5 
NG = NG + I 
IFING.LE. 10) GO TO 4 

5 CALL RESUL i ( NNOD ,V ,EI , lER, ITMAX, NCODE ) 

STOP 

END 
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SUBROUT INE STORE ( NDOF , NBCt TL ,NEL EM , 

1 KATEST,K£TEST,kQTEST ,KA,KEI , 

2 KQ,KS,EPS,MSIG, 

3 ASAV, WSAV,EISAV,P0,SS,SK, ICCRR) 

FOR LINEAR PROBL EMS , OM I T 
•IMPLICIT REAL^i'S ( A-H, 0-Z) • CARD 

IMPLICIT Rt AL«8( A-H,0-Z) 

DIMENSION ASAV(IO) tCSAVdO) ,£ISAV(10),SS(10), 

1 SKI 10) , Pj( 10) , ICORRI 20,4) 

THE FOLLOWING PARAMETERS ARE USED FOR 
TESTING THE CONVERGENCE OF THE PROGRAM. 

NELMAX=8 

KETEST=l 

KQTEST=1 

KATEST=1 

EPS=1.0D-06 

NSIG=5 

READ IN ALL INFORMATION NEEDED TO THE 
ENTIRE FIELD OF THE PROBLEM. 

READ(5,lu) KA,KEI ,KQ,NELEM 
REA0(5,3) NDOF , NBC ,KPO,KS ,TL 

3 FOPMAT(4I5,G15.5) 

READ(5,4) (POm ,I = l,KPO) 

READ(5,4) (ASAVI I ) ,1=1 ,KA) 

REAOI5,4> (QSAVI I ) , i = l ,KQ) 

READ(5,4) (EISAVI I ) , 1=1 ,KEI ) 

READ(5,4) (SKI I ) , I =1 ,KS) 

READ(5,4) (SS( I ) . 1=1 ,KS) 

4 F0RMAT(4G15.3) 

READ IN THE CORRESPONDENCE TABLE 
BETWEEN LOCAL AND GLOBAL COORDINATES 

DO 11 I=1,NELMAX 

11 REAO(5, 10) ICORRI 1,1), 1C0RR( 1 , 2 ) , ICORRd , 3) , ICORRI 1,4) 
10 F0RMAT(4I5) 

RETURN 

END 
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SUBROUTINE I NCHK ( NEL EM , NNOD , NBC.NOOF, 

XtTLiEI ,Q,IP, JP,KP,ISOIL, 
SSOlLtALPHA,NCODE tVSCALF, 

ASAV, EISAV.QSAV, SK.,SS, IDBC , V,PO) 

FOR LirJEAR PROBLEMS, OMIT 
•IMPLICIT REALMS ( A-ri, O-Z) • CARD 

IMPLICIT REAL-0( A-H,0-Z) 

DIMENSION ASAV(IO) tEISAVI 10) ,QSAV( 10) , 

1 IDBCIZO) ,V(20),PO(10),SK(10),SS(10) 
VSCALE= 1. 0-06 
Q=QSAV( IP ) 

•EI=EI SAVUP) 

IFINCODE.NE. 1 ) GO TO 7 

THESE TWO FOLLOWING STATEMENTS IS DEALING 
WITH NATURAL SOILS 

SSCIL=SS( ISOIL) 

ALPHA=SKI ISOIL) 

GO TO 9 

SSCIL MAY HAVE A CERTAIN VALUE IF THE 
SHEAR FORCE IN THE FOUNDATION MATERIAL 
IS TAKEN INTO ACCOUNT. 

7 SS0iL=0.0 
ALPHA=ASAV(KP) 

9 NN0D=2* (NELEM+1 ) 

X=TL/{ 2*NELEM) 

PRINT 6,NNUD,NELEM,NBC,N00F,X 

8 FORMATISX,* NUMBER OF NODAL PO I NTS = ‘ , I 5 , / , / , 
15X, ‘NUMBER OF EL EMENT S= • , I 5 , / , / , 

2 5X, 'NUMBER OF BOUNDARY COND I T I ONS= ' , I 5 , / , / , 
35X, ‘NUMBER OF DEGREES OF FREEDOM= • , I 5 , / , / , 
45X, • ELEMENT LENGTH-’, IX, 015. 6,////) 

B£TA={ ALPHA/ (4. -El ) )*=«=. 25 

THE FOLLOWING BOUNDARY CONDITIONS MAY NOT 
BE NECCESARY IN GENERAL AND CAN BE 
READJUSTED IN THE SUBROUTINE BOUND. 

IDBC(1)=1 

I CBC( 2 ) =NNOD 

V(1)=0.0 

V( NNOD) =0.0 

RETURN 

END 
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SUBROUTINE STIFF2 
IMPLICIT REAL*8( A-H,0-Z) 

C 

C- THIS SUBROUTINE IS USED TO COMPUTE 
C THE ELEMENT STIFFNESS MATRICES AND THE 
C SYSTEM STIFFNESS MATRIX FOR NONLINEAR 

C PROBLEM , ACCORDING TO THE FORMULA ( ) 

C ON PAGE 

C 

COMMON EK(4,4) ,DC(4,^),0M(4t4) ,DN(4,4,<V) , 

1 IC0RR(20,4) , SYSK( 20, 20 ) ,SYSC( 20,20t ZOI , 

2 SYSF(20),EF(4), IDBC(20),B(10) , 

3 NNOD, NELEM, NDOF ,NBC , T L , AA 1 , AA2 , 

4 X, El, Q, ALPHA, SS01L,Al,A2 
DO 5 I=1,NNC0 

DO 5 J = I , N NOD 
SYSK( I , J) =0.00+0 
DO 5 K=1,NNUD 
5 SYSC( I , J,KJ=0.0D+0 
DO 6 N=1 , NELEM 
C 

C THE FOLLOWING FUNCTIONS SFUNCT , EFUNCT , 

C B2 , B3 ARE DEALING WITH VARIABLE SHEAR 

C FORCE IN THE FOUNDATION MATE RI AL S , VARI AB LE 

C BEAM FLEXURAL RIGIDITY AMD NONLINEAR 

C REACTION OF THE FOUNDATION . THE REACTION 

C OF THE FOUNDATION MAY BE WRITTEN AS 

C R (X) = 62- V + 33-V**2 

C 

SS01L=SFUNCTIN,X) 

EI=EFUNCT(N,X) 

B( 1)=0. OD + 0 
B(2)=S2 (N, Xj 
B(3)=B3( N, X) 

SSOIL=SSOIL/EI 
A1=ALPHA/EI 
AA1=A1*B( 2) 

C 

C THE FOLLOWING IS THE BEAM ELEMENT 

C STIFFNESS MATRIX 

C 

EK( I, 1 ) =12.D+00*( I .0+00/X^*3) 

EK(i,2)=O.D + 00*( I .D+00/X=i==!=2 ) 

EK( I ,3) =-12.D+00-( l.D+00/X**3) 

EK( l,4}=6.D + 00^( 1.0 + 00/X>!'«2) 

EK(2,2) =^.0 + 00=!^( 1 .D + OO/X) 

EK(2,3) =-6.0 +00* ( 1.0+00/X**2J 
EK(2,4)=2.0+00v( I .D+OO/X) 

EK(3, 3) =12 .D+00* ( 1 .D+OO/X** 3) 

EK(3,4) =-6.0+00* ( 1 .D+OO/X** 2) 

EK( 4, 4 1 =4.0 + 00*1 1 .0 + 00 /X) 

EK(2,n=EK(l,2) 

EK(3,1)=EK(1,5) 

EK(3,2)= £K(2,3) 

EK(4, 1) =EK( 1 ,4) 

EK(4,2)=EK(2,4) 

EK(4,3)=EKI3,4) 
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THE FOLLOWING IS THE ELEMENT FOUNDATION 
STIFFNESS MATRICES , ASSOCIATED TO THE 
LINEAR TERMS 

0C(1 , li =1 .20<-00/X 
DC! 1,2) =0.10+00 
0C(l,3)=-0C(l,l) 

DC( 1 ,4) =DCU ,2 ) 

DC(2,1 )=DC(1 ,2) 

DC(2,2)=0.133333D+00*X 
DC(2,3) =-DC( 1,2) 

DC(2,4)=-0.033333D+00«X 
DC(3,1) =DC(1 ,3) 

DCI3,2)=DC(2,3) 

OC(3,3)=OC(l,l) 

0C(3,4)=DC (2,3 ) 

DC(4,1)=DC(1,4) 

DC(4,2>=0C (2,4) 

DC(4,3)=DC(3,4) 

DC(4,4)=DC(2,2) 

DM (1,1 ) =.371429D+00*X 
DM(l,2) = .052381D+O0^^X=!‘=^2 
DM( 1 ,3 )=.123571D+00*X 
DM (1,4) =-. 030952 0 + 004X=:=*2 
DM! 2,2 ) = .0095240 + 0 0-X=!‘ =4=3 
DM( 2,3) = .03095 2D + 00«X*=!'2 
DM( 2,4) =-.007143D+00*X*«3 
DM(3,3) = .37i4290+00='=X 
DM(3,4) =-.052381D+00*X**2 
0M(4,4) = .009 5240^ 00=i=X«=)'3 
DM(2, 1 ) =DM( I ,2) 

DM(3,1) =DM(1,3) 

DM(4, 1 ) =DM( 1,4) 

DM(3,2 )=DM(2,3) 

0M(4,2) =DM(2,4) 

DM(4,3) =DM(3,4) 

THE CORRESPONDING COMPONENTS OF ALL 
STIFFNESS MATRICES ASSOCIATED TO THE LINEAR 
TERMS MAY BE ADDED AS FOLLOWS 

DO 1 I=1,ND0F 
DO 1 J=1,ND0F 
0C( I , J) =SSOIL*DC ( I , J) 

DM(I , J) =AA l=i=DM( I ,J) 

1 EK(I,J)=EK(I,J)+OC(I,J)+DM< I,J) 
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THE FOLLOWING ARE THE COMPONENTS OF THE 
MATRIX ASSGCIATcO TO THE NONLINEAR TERMS 

ON ( 1, 1, 1 ) = .30 7143D4-00«X 
ON(2tl, 1 ) = .03 849?D+-00*X*#2 
ON (3, 1, 1 )=. 0642 360 +00*X 
0N(4t 1 , U = -.017061D+00*X=i'*2 
ON(2f2f I ) = . 0065490 fOO=!‘X=:'#3 
ON (3,2, 1 )=.013 839D4-00>i^X**2 
ONI 4 1 2,1 )=-.0035 71D + 00=i‘X*=!'3 
0N(3,3, l)=.0642860fOO*X 
ON (4, 3, 1J=-.OI3889D+00*X**2 
ON (4, 4, 1 ) = .003l 750 + 00^-X=!'*3 
ONI 2,2, 2) = .00 11910 + 00*X^=!'4 
ONI 3,2, 21= .0031 750 + 00* X**3 
0NI4,2, 2 )=-.0007940+00*X**4 
0NI3,3, 2)=.0170610+00*X**2 
0NI4,3, 2 J=-.0035 71D + 00*X**3 
0NI4,4, 2)=.0007940+00*X**4 
DNI3,3,3 )=.:)071430 + 00*X 
ONI 4,3, 3 ) = -.038492 0 + 00*X**2 
ONI4,4,3)=.0063490+00*X**3 
0NI4,4, 4)=-.0011 910+00*X**4 
ONI1,2,1)=OMI2,1,1) 

ONI 1,3, 1 )=0M(3, I , 1 ) 

0NI2,3,1)=DN{3,2,1) 

0NIl,4,l)=DN(4,l,n 
0NI2,4, 1 )=0NI4, 2, 1 ) 

0NI3,4,1 )=DM(4,3 , 1 ) 

ONIl,l, 2 )=DNI 2 ,l,n 

ONI 1,2, 2) =ON I 2,2,1 ) 

ONI 1, 3, 2J=ON(3,2,l ) 

ONI l,4,2)=ON(4,2, 1 ) 

0NI2,1,2)=DN(2,2,1 

0NI2,3,2)=DN(3,2,2) 

ONI2,4,2)=ON(4,2,2J 
0NI3,1,2»=0NI3,2,1 ) 

0NI3,4, 2)=0NI4,3 ,2 ) 

0NI4, 1 ,2)=0N(4, 2, 1 ) 

ONI 1,1, 3 ) = ON 1 3,1 ,1) 

ONIi,2,3)=ON(5,2,n 
ONI 1,3,5 )=0N( 3,3, U 
DNI1,4,3) = 0NI4,3,1 ) 

ONI2, 1,3)=0NI3,2,1 J 
ONI2,2,3)=OM(3,2,2) 

0NI2,3,3)=0N(3,3,2) 

DNI2,4,3J=0N(4,3,2) 

ONI3,l ,3) = 0N(3,3,1 ) 

ONI3,2,3J=ON(3,3,2) 

ONI 3,4 , 3 ) = 0NI 3 ,3 ) 

0NI4, 1, 3)=DN( 4,3, 1) 

0NI4,2,3)=0NI4,3,2) 

ONIl, 1,4)=0N(4,1,1 ) 

ONI l,2,4)=0N(4,2,U 
0NIl,3,4J=DNI4,3,l) 

ONI 1 ,4,4) = 0N(4,4, 1 ) 

ONI2, 1,4) = 0NI4,2,1 ) 

0NI2,2,4)=0NI4,2,2) 

0NI2,3,‘t)=DN(4,3,2) 

ON 12, 4,4)= ON (4, 4, 2) 

0NI3, 1, 4)=0N( 4,3 ,1 ) 

ONI3, 2,4) = QNI4,3,2 ) 

ONI3,3,4)=ON(4,3,3) 

0NI3,4,4)=0NI4,4,3) 

DN(4, I, 4) = DN(4,4, 1 ) 

0N(4,2,4)=0N(4,4,2) 

0NI4,3,4)=DNI4,4,3) 



89 



c 

c 

c 

c 

c 



THE SYSTEM STIFFNESS MATRIX IS OBTAINED 
BY COLLECTING THE ELEMENT STIFFNESS 
MATRICES OVER ALL ELEMENTS OF THE SYSTEM 



AA2=Ain«B(3) 

DO 4 K=l,NDOF 
DO 4 J=l,NOOF 
DO 4 I=l,NDOF 
4 DNU t Jt K)=Aa 2*DN(I , JtK) 

DO 7 I=1tNDOF 
NSl= iCORR(N,n 
DO 7 J=ltNOOF 
NS2= ICORRINtJ) 

SYSK(NS1,NS2) = SYSK(NS1,NS2) + EKdfJi 
DO 7 K=1 tNDOF 
NS5= ICORR(N,K) 

7 SYSCI NSl,N52,NS3)=SYSC(NSL,NS2tNS3) + DN(l-,J,K) 
6 CONTINUE 
RETURN 
END 
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SUBROUTINE LOAD 



FOR THE LINEAR PROBLEMS , USING SINGLE 
PRECISION, THE 'IMPLICIT RE AL-8 ( A-H,0-Z ) * 

CARD MUST BE OMITTED 

IMPLICIT REAL*8( A-H,0-Z) 

THIS SUBROUTINE IS USED TO COMPUTE 
THE ELEMENT FORCE VECTOR 
AND THE SYSTEM FORCE VECTOR 

COMMON EK( 4,4) ,DCI 4,4) ,DM(4,4J ,DN(4,4,4) , 

1 ICCRR(20,4) ,SYSK(20,20),SYSC(20,20,20), 

2 SYSFI 20 ) , EF (4 ) , ID3C(20) ,8( 10) , 

3 NNOD,NELEM, NDCF ,NBC , TL , AAl , AA2, 

4 X, El, Q, ALPHA, SS0IL,A1,A2 
DO 1 I=l,NNOD 

1 SYSFI 1 )=0, 00+00 
DO 3 N=l,NELEM 
EF(1)= .5D+00«X 
EF<2)= ,083333D + 00=s=X**2 
EFI3)= ,5D+00^=X 

EF(4)= -.033333D+00v-X*«2 
Q=QFUNCT(N,X)-A1-8(1 ) 

EF(U=Q-'‘EFa) 

EF(2)=Q=^EFI2) 

EF(3)=0-£F(3) 

EF(4)=Q«EF (4) 

DO 2 NT=1,ND0F 

2 EF(NT)=EF(NT)/EFUNCT(N,X) 

I=ICORR(N, 1) 

J=ICORR(N, 2) 

K=ICORR (N, 3) 

L=IC0RR(K,4) 

SYSF I I ) =SYSF( I ) + EF( 1 ) 

SVSFIJ) =SYSF(J) + EFI2) 

SYSFI K) =SYSF(K) + EFI3) 

3 SYSFI L)=SYSFIL ) + EFI4) 

RETURN 

END 
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SUBROUTINE SOUND 



fOR LINEAR PROBLEMS, THE 'IMPLICIT REALM'S 
(A-H,0-Z)' CARD MUST BE OMITTED 

IMPLICIT REAL<=8( A-H,0-Z) 

BEFORE SOLVING THE ALGEBRAIC EQUATION SYSTEM 
THE BOUNDARY CONDITIONS MUST BE APPLIED 

COMMON EK(4,4) ,DC( A,4) ,DM(A,4) ,DN(4,4,4) , 

1 I CORK (2 0,4) ,SYSK( 20 , 20 ) , S YSC ( 20 , 20 , 20 ) , 

2 SYSFt 20) , EF(4) , I06C(20) ,B( 10) , 

3 NNODtNELEM, NDOr ,NBC , TL , AAl , AA2, 

4 X, El, Q, ALPHA, SS0IL,Al,A2 
DO 3 N=1,N8C 

FOR FREE-FREE ENDS BEAM IT IS DESIGNED 
TO SKIP THE FIRST BOUNDARY CONDITION 
WHICH IS RELATED TO THE DEFLECTION AT THE 
END . THE SLOPE AT THE MIDDLE , WHICH 
IS ZERO BY SYMMETRY, MUST BE IMPOSED. 

FOR OTHER TYPES OF BOUNDARY 
CONDITIONS, THIS SUBROUTINE IS NOT TRUE. 

IF<-N.EQ,U GO TO 3 
KK=IDBC (N) 

DO 4 K=l,MNOD 
SYSK( KK,K) =O.CD<-00 
SYSK{K,KK)=0.00 4-00 
DO 4 I=1,NN0D 
DO 4 J=1,NN0D 
SYSC( KK, I , J) =0.004-00 
4 SYSCIKK , J, I ) =0.00+00 
SYSC( KK,KK,KK) =0.00+00 
SYSF(K,K)=0.0D + 00 
SYSK<KK,KK)=l.D+00 
3 CONTINUE 
RETURN 
END 
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SUBROUTINE CURE I T { LP , PO, B ) 

USING LSQPL2 LIBRARY SUBROUTINE , THE 
IRRATIONAL POWER CURVE WITH P NOT EQUAL TO 
AN INTEGER IS REPLACED BY A LEAST SQUARE 
BEST FIT SECOND ORDER POLYNOMIAL 

IMPLICIT REAL*8( A-H,0-Z) 

DIMENSION 6( 10) ,WI (100) ,Y{ 100) ,DELY( 100) , 

1 P0( 10) ,SB( 100) , F2( 100) ,X( 100) ,ER( 100) 

REAL’i^S TITLE(10)/10-' '/ 

DATA WI/100*1.D0/ 

NMAX=-2 
P=PO(LP ) 

IF(P.EQ.2.D+00) GO TO 10 
N=5 

DO 2 1=1, N 
FLOATI=I 
X( I ) = FLCATI 
2 F2( I )=X( I ) v*p 

CALL LSGPL2(N,NMAX,X,F2,WI,Y ,DELY,B,SB, TITLE) 

DO 6 1=1, N 

IF(F2( I ).LE.l.D-06) GO-TO 7 
ER( I )=D6LY( I )/F2( I ) 

GO TO 6 

7 ER(I)=0.0D+00 
6 CONTINUE 

ERR0R=ER(1) 

DO 8 1=2, N 

IF(OABS(ER( I )) .LE.ERROR) GO TO 8 
ERR0R=DA6S(ER{ I) ) 

8 CONTINUE 
PRINT 9, ERROR 

9 FORMAT!//, lOX, • THE MAX RELATIVE ERROR IS=',D15.7) 
GO TO 11 

10 B( l)=0.0D«-00 
B(2)=0.0D<-00 
B(3)=l .D+00 

11 RETURN 
END 
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SUBROUTINE FSOILIB) 

I^PLICIT REALMS! A-H,0-Z) 

THE EXPERIMENTAL DATA ON A TYPICAL 
NATURAL SOIL, I. E . .REACTION OF FOUNDATION 

VERSUS LOADING CONDITION FROM TABLE 

AND FIGURE ,ARE APPRIMATED BY A LEAST 

SQUARE BEST FIT SECOND ORDER POLYNOMIAL 
USING LSQPL2 LIBRARY SUBROUTINE 

DIMENSION B( 10) , WI ( 10) ,Y( 10) ,DELY( 10) , 

1 SB( 10),E2( 10) ,X( 10) ,ER(10) 

REALMS TITLE! 10)/10=^' •/ 

DATA WI/10=i'l.D0/ 

NMAX=-2 

N=6 

DO 2 1 = 1, N 

2 READ! 5,3) X( I ) ,F2( I ) 

3 FORMAT! 2G15.6) 

CALL LSQPL2!N,NMAX,X,F2,WI , Y,DELY,B, SB.TITLE) 

THE SEA SEDIMENT FOUNDATION IS MODELED 
FROM INLAND SOIL DATA BY A SCALE 1:3000. 

OMIT THE Three following cards in case 
IN CASE OF INLAND FOUNDATION 

B ! 1) = 6! 1 )*8.3D+-0 0/30 00.D+00 
B( 2) = 8! 2)=i'8.3Df 00/3000.Df00 
B!3) = B! 3) 3D ♦■00/3000, D+-00 

RETURN 
END 
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SUBROUTINE VGUES S ( NNOD, NELEM, VSCALE, ITMAX,VJ 

THIS SUBROUTINE PROVIDES AN INITIAL 
ESTIMATE FOR NONLINEAR PROBLEMS 

IMPLICIT REAL*3( A-H,0-Z) 

DIMENSION V(20) 

ITMAX^lOO 

THE MAIN PROGRAM IS DESIGNED IN 
SUCH A WAY THAT IF ONE INITIAL ESTIMATE 
DOES NOT GIVE SOLUTIONtTHE SCALE WILL BE 
CHANGED, AND THE ITERATION RESTART AGAIN. 

HOWEVER , AFTER 10 CYCLES, IF MO SOLUTION IS 
OBTAINED, THE READER MAY CHANGE THE INITIAL 
SCALE BY HAND IN THE SUBROUTINE INCHK 

VSCALE=1 . 5D+00*VSCALE 
TAU=0. OD+00 
NT=NELEM+l 
DO 22 1=1, NT 
11 = 2*1 
111 = 2 * 1-1 

FOR SIMPLICITY , A SECOND ORDER POLIMOMIAL 
IS USED FOR ESTIMATE SOLUTION . THE READER 
MAY DESIGN OTHER ESTIMATE FUNCTION 
AS HE WHISHES 

V( I I ) = VSCALE*2.D<-00*(4.D<-00-TAU) 

V ( II I ) =VSCaLE*( ( 8 .D + 00*TAU-TAU**2I <-G.01D*-00) 

22 TAU=TaU<-1 .D+00 

vm=o.oD+oo 

PRINT 22, (V( I ) , I=l ,NNOD) 

23 FORMAT! ////,15X, • INITIAL GUESS :’,///,( 1 5X , D1 5 .6 ,/ )) 
RETURN 

END 
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SUBROUTINE DIFSOL 



USING FINITE DIFFERENCE METHOD TO SOLVE 
NONLINEAR PROBLEMS , Tri I S SUBROUTINE IS 
USED STRICTLY FOR THE CASE OF UNIFORM BEAMS 
WITH UNIFORM LOADS AND THE REACTION OF 
FOUNDATION IS PROPORTIONAL TO THE SQUARE 
OF ITS DEFLECTION 

IMPLICIT REAL*8( A-H,0-ZJ 

COMMON EK(4,4) ,DC(4,4) ,DM( A,A) ,ON(4,A,A) , 

1 IC0RR(20,4) . SYSK( 23,20 J tSYSC I 20,20, 20) , 

2 SYSFI 20) ,EF(4) , I0BCI20) ,B( 10) , 

3 NN0D,NELEM,ND0F,NBC,TL,AAl,AA2, 

4 X, El, Q, ALPHA, SS0IL,A1,A2 
DIMENSION U( 10) ,W0( 100) 

TWO BEAM MODELS ARE CONSIDERED : 

EXTERNAL DIFF4 CORRESPONDS TO 4-ELEMEMT 
MODEL, AND EXTERNAL DIFF8 TO 8-ELEMENT MODEL 

EXTERNAL DIFF4 
EXTERNAL DIFF8 
I FORMAT! • 1* ) 

ND=NELEM 
OD=ND 
NSIG=5 
EPS=l.D-Oe 
NCCUNT= 1 

USCALE=1000.D+00 

3 ITMAX=100 
NC0UNT=NC0UNT+1 
TAU=0. OD+00 
USCAL£=1. 5Di-00=<-USCALE 
DO 4 I=1,ND 
TAU=TAU+2 .DfOO/DD 

4 U! I)=U5CAL£*(TAU-. 25D+00*TAU««2) 

IF(ND.GT.2) GO TO 5 

CALL ZSYSTM! DIFF4 , E PS , NSI G , NO , U , ITMAX, WD ,PAR, lER) 
GO TO 8 

5 CALL ZSYSTM(DIFF8, EPS, NSIG,ND,U, ITMAX, WD, PAR, lER) 

8 IF(NCOUNT.GT.30) GO TO 6 
IF(IER.NE.O) GO TO 3 

6 PRINT 7 , lER, ITMAX 

7 FORMAT! ////, 15X, • I ER= • , I 5//15X , • ITMAX=* ,15) 

DO 9 1=1, NO 

9 U! I )=U! D/EI 
PRINT 10 

10 FORMAT! ////, 15X, 'F.D.M. RESULT',////) 

PRINT 1 1,!U! I) ,1=1 ,ND) 

11 FORMAT ( 15X,G20. 6/) 

RETURN 

END 
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SUBROUTINE STIFFl 

C 

C THIS SUBROUTINE IS USED TO COMPUTE THE ELEMENT 

C STIFFNESS MATRICES AND THE SYSTEM STIFFNESS 

C MATRIX FOR LINEAR PROBLEMS 

C 

CCMHON EK(4,4),DC(4,4),DM('V,4)iDN(A,4,4J, 

1 IC0RR(20, 4) ,SYSK( 20,20) tSYSC(20, 20,20) t 

2 SYSFI 20) , EF(4) , IDBCI 20) ,B( 10) , 

3 NNOD, NELEM, NDOF ,NBC , TL , AA 1 , A A2 , 

4 X, El, Q, ALPHA, SS0IL,A1,A2 
ALPHA=ALPHA/EI 
SSCIL=SSOIL/EI 

EKI i, 1) = 12.*(1./X*«3) 

EK( 1,2)=6.*( l./X-*2) 

EK(1,3)=-12.*(1./X=«=*3) 

EKI 1 ,4)=6.=<=( i./X<^=i‘2) 

EK(2,2)=4.*( l./X) 

EK(2,3)=^6.*( l./X**2) 

EK(2,4) = 2.v-( l./X) 

EK(3,5)=12.-(1 ./X*=i'3) 

EK(3,4) =-6.«( 1 ./X=!'=i‘2) 

EKI4, 4) =4.’!=( l./X) 

EK(2,1 )=EK(1,2) 

EK(3,l)=EKa,3) 

EK(3,2)= £K(2,3) 

EK(4,1)=EK(1,4) 

EK(4,2)=EKt2,4) 

EK(4,3)=EK(3,4) 

C 

C THE FOLLOWING ARE THE ELEMENT STIFFNESS 

C MATRICES ASSOCIATED TO THE FOUNDATION 

C 

DC( 1 , 1 ) =1 .2/X 
DC (1,2) =0.1 
0C(l,3)=-DC(l,l) 

0C(1,4)=DC(1,2) 

DC(2,1) =DC( 1,2) 

DC(2,2) =0.133333*X 
DC(2,3) =-DC( 1,2) 

DC(2,4)=-0.G333:k3-'!'X 
DCI3, 1) =DC( 1 ,3 ) 

DC(3,2)=DC(2,3) 

0C(3,3)=DC(1,1 ) 

DC(3,4)=DC(2,3) 

0CI4, 1) =DC( 1 ,4) 

0C(4,2)=DC(2,4) 

DC(4,3) =DC(3,4) 

DC(4,4) =DC(2,2) 

DM( 1, 1 )=.371429«X 
OM(l,2)=.05238l«X*«2 
DM(l,3)=.12d571*X 
0M(l,4)= -.030952*X=<'*2 
DM(2,2 ) = .009524=f^X*«3 
DM(2,3)= .030952*X«’?2 
DM(2,4) =-.007143*X=!‘*3 
DMI3,3) = .371R29*X 
DM(3,4)= -.052381=!'X**2 
DM(4,4)= .G09524*X**3 
DMI2 , 1 ) =0M( 1,2) 

0M(3, 1)= DM! 1,3) 

DM(4, 1 ) = DM( 1, 4) 

OM(3,2)= DM(2,3) 

DM(4,2)=DM(2,4) 

0M(4,3)= DM(3,4) 
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SINCE ALL ELEMENT STIFFNESS MATRICES ARE 
ASSOCIATED TO THE LINEAR TERMS t THEY 
MAY BE ADDED AS FOLLOWS 

DO 1 1= l,NDOF 

DO 1 J= l.NDOF 

1 EK(I,JJ=EK( I , JJ+ALPHA^i^DMI I , J ) «-SSOI L*DC (I ,J) 
DO 2 I=1,NNGD 

DO 2 J =l,NNOD 

2 SYSKI I , 0.0 

THE SYSTEM STIFFNESS MATRIX IS OBTAINED BY 
COLLECTING THE STIFFNESS MATRICES OVER ALL 
ELEMENTS OF THE SYSTEM 

DO 3 N= UNELEM 
DC 4 i=l,NDOF 
NS1= ICORR(N,I) 

DO 4 J= 1 , NDOF 
NS2= ICCRR(N,J) 

4 SYSKI NSi,NS2J= SYSKI NSl ,NS2 ) + EKII,J) 

3 CONTINUE 
RETURN 
END 
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SUBROUTINE SOLVE ( SYSK , SYSF , V , NNOD ) 

BASED ON GAUSSIAN ELIMINATION TECHNIQUE 
THIS EQUATION SOLVER IS USED TO SOLVE 
THE LINEAR ALGEBRAIC SYSTEM ONLY 

DIMENSION SYSF(20) ,SYSK(20,20) , 

I IV(2D) ,SYSFP(20) ,V(20I 

DC 1 K=l,NNOD 
1 IV(K)=K 



LlOO 



1300 

1200 



1006 

1005 

1004 



11 

10 



GAUSSIAN ELIMINATION 

NMl=NNOD-l 

DC 1004 1=1, NMl 

IP1=I+1 

ATGP=A6S(SYSK( 1,1)) 

IPRIME=I 

DO 1100 IP=IP1,NN0D 

IF(ABS( SYSK( I?,l ) ) .LE.ATCPJ GO TO 1100 
AT0P=A8SI SYSK( IP , 1 ) ) 

I PRIHE= I P 
CONTINUE 

IF( IPR IME.EO.I ) GO TO 1200 
ISAV=IV ( IT'RIME) 

IV( IPRIME)=IV( 1 ) 

IV( I )=I SAV 
DO 1300 J=l,NNOD 
TEMP=SYSK( IPRIME , J ) 

SYSK( IPRIME, J) =SYSK( I , J) 

SYSK( I , J1=TEMP 
DO 1005 K=IP1,NN0D 
FACTOR=SYSK(K, I ) /SYSK( 1,1) 

SYSKIK, I ) =FACTQR 
DO 1006 J=IP1,NN0D 

SYSKIK, J)= SYSK (K,J)-F ACTOR=^= SYS KI I , J) 

CONTINUE 

CCNT INUE 

CONT INUE 

DO 9 1=1, NNOD 

K=IV( I ) 

SYSFP( I )=SYSF( K) 

DC 10 J=1,NM1 
JPl=Jf 1 

DO li K=JP1,NN00 

SYSFPI K)=SYSFP{K )-SYSK( K, J )*SYSFP( J ) 

CONTINUE 

CONTINUE 



BACK SOLUTION 

SY SF( NNOD) =SYSFP( NNOD) /SYSK (NNOD, NNOD l 

DO 13 I 1 = 1) NMl 

I =NNOD-I I 

IPl=l*l 

SUM=0.0 

DC 14 J=IP1,NN00 
14 SLM=5UM+SYSK( I , J )#SYSF( J) 

13 VU ) = ( SYSFPI 1 )-SUM)/SYSK( I , I) 

RETURN 

END 
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SUBROUTINE RESULT ( NNOD , SY SF , El , I ER, ITHAX , NCODE ) 



THIS SUBROUTINE IS USED TO PRINT THE 
SOLUTION UF EITHER LINEAR OR NONLINEAR 
PROBLEMS. FOR LINEAR PROBLEMS , OMIT 
•IMPLICIT REAL<^8( A-H, C-Z) • CARD 

IMPLICIT REAL«8( A-H,0-Z) 

DIMENSION SYSF(20) ,V( 20) 

IF (NCODE. EQ.IJ GO TO 19 
PRINT 28 , IER,ITMAX 
28 F0RMAT(//,15X,» I E R = » , IX , I 2 , /// / , 
115X,MTMAX =• , IX ,I3f////) 

19 NN=NN00/2 

DC 20 I=1,NN 
I 1 = 2*=I 
JJ=2=«^i-l 

20 V( IJ=SYSF( JJ) 

PRINT 30 

30 FORMAT! ////t I5X, • DEFLECTION :*,//) 

PRINT 31, ( vn ) , I =1 ,NN) 

31 FORMAT! 5X,G20.6,/) 

IF THE DEFLECTIONS AS WELL AS THE SLOPES 
AT THE NODAL POINTS ARE DESIRED , 

ADD THE FOLLOWING CARO 
PRINT 31, !SYSF< I ) , I = 1,NN0D) 

RETURN 

END 
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THE FUNCTION AUX WHICH IS TRANSLATED 

FROM EQUATION PAGE .GENERALLY 

REMAINS THE SAME. BUT THE READER HAS TO 
REDESIGN OTHER EXTERNALS , SUCH AS QFUNCT, 
EFUNCT,SFUNCT,BZ,B3,WHICH MAY VARY , 
DEPENDING CN LOADING CONO I T ION S . BE AM 
CHARACTERISTICS AND TYPES OF FOUNDATIONS 



DOUBLE PRECISION FUNCTION AUXIV.K.PAR) 

IMPLICIT REAL*8( A-H.O-Z) 

CCMMGN EKIA.4) ,DC(4,A),DM(A,4) ,DN(At4,4) , 

IC0kR(20,4) , SY SKI 20,20 ) .SYSCI 20,20, 20 J , 
SYSFI 20),EF (4), I0BC(20) ,6(10), 

NNOD.NELEM, NDOF ,N8C , TL , AAl , AA2, 

X, EI,0,ALPHA,SS0IL,A1 ,A2 
DIMENSION V(20) 

SKTEMP = 0 .004-00 
SCT£MP=0.0D+00 
DO 1 i=l,NNCD 

SKTEMP=$KTEMP4-SYSK(K,I )*V( I) 

DO 1 J=I,NNOD 

1 SCT£MP=SCTEMP+SYSC(K, I , J ) « V ( I ) =i=V ( J ) 
AUX=StCTEMP+SCTEMP-SYSF(K) 

RETURN 

END 

DOUBLE PRECISION FUNCTION QFUNCT(N.X) 

IMPLICIT REAL*8( A-H,0-Z) 

FLOATN=N 

XI=(FL0ATN-.5D+00)*X 
IF(N.NE.l) GO TO 3 
QFUNCT=0. 00+00 
GO TO 4 

3 QFUi\CT=((3.5D + 00/3.0+00) + ( 1 .0+0 0/ 150.D + 00) * 

1 XI)’!^l.p+04 

4 RETURN 
END 

DOUBLE PRECISION FUNCTION B2(N,X) 

IMPLICIT REAL#3( A-H,0~Z) 

FLOATN-N 

XI=( FLCATN-.5D+00i«X 
IF(N.NE.l) GO TO 3 
B2=2000.0+00+40.D+00#XI 
GO TO 4 

3 82=( 100.D+0U/15.D+00) «XI + (11000.0+00/3.0+00) 

4 RETURN 
END 

DOUBLE PRECISION FUNCTION B3(N,X) 

IMPLICIT REAL-31 A-H,0-Z) 

FLOATN=N 

XI = (FL0ATN-.5D + 00)=4‘X 
B3-.5D+00«Xi +300 .0+00 
RETURN 
END 

DOUBLE PRECISION FUNCTION EFUNCT(N.X) 

IMPLICIT REAL*8( A-H,0-Z) 

FL0ATN=N 

X 1={ FL0ATN-.5D+00)#X 
IF(N.Nc.l) GO TO 3 
EFUNCT=1 .D+ll+( 1 .0+11/50.D+00) *XI 
GO TO 4 

3 EFUNCT=2.D+11 

4 RETURN 
END 

DOUBLE PRECISION FUNCTION SFUNCT(N,X) 

IMPLICIT REAL*8( A-H,0-Z) 

SFUNCT=0.0D+00 

RETURN 

END 
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DOUBLE PRECISION FUNCTION DIFF4{ U, K, PAR) 

IMPLICIT REALMS! A-H,0-Z) 

COMMON EK(4,4) ,DC(4,4) »0M(4,4) ,0N(4,4f4) , 

1 IC0RR<Z0,4) ,SYSK( 20t20) ,SYSC( 20t20t20) , 

2 SYSFI 20) ,EF (4 ) ,I DBC(20) ,B( 10 ) , 

3 NNCD, N5LEM, NDGF ,MSC, TL, AAl » AA2, 

4 X»EI ,Q, ALPHA, SSOILtAl ,A2 
DIMENSION U(10) 

GO TO (5,10) ,K 

5 DIFP4=1 . 536D+OC*U(l)+0.001D + 00*A2=i'(TL=i=*4)=^‘U(l)**2- 
1 1.024D + 00*U( 2)-0.001D + 00=i-'Q=!=(TL=<=«4) 

RETURN 

10 DIFF4=-2.048 0f00=«'U( 1 ) + 1.536D+00*U( 2) + 

10. 00ID^-00^A2*( TL**4) -U( 2)*=^2-0. 00 i0*00« Q* ( TL**4J 
RETURN 
END 



DOUBLE PRECISION FUNCTION D I FF 8 ( U, K, PAR ) 

IMPLICIT R£AL=!=3( A-H,0-Z) 

CCMMGN EK(4,4) , DC (4,4) ,DM( 4,4) ,DN( 4,4,4) , 

1 IC0RR(20,4) , SYSK(20, 20) ,SYSC(20, 20,20) , 

2 SY5F(20) ,£F{4) ,IDBC(20) ,B( 10) , 

3 NNOD, NELEM, NDOF ,NBC , TL , AA 1 , AA2 , 

4 X, El, Q, ALPHA, SS0IL,A1,A2 
DIMENSION U(10) 

GO TO ( 5,10,15,20) ,K 

5 0IFF8=5.0+G0*U( 1 )-4.D+00*U( 2) *■ U(3) + 

1 (X=!^=!-'4) 4-( A1#B( l)+AAi-U( 1) +AA2-UI 1)*42) 

2 -q*x**4 
RETURN 

10 DIFF8=-4 ,D*-00-U( 1 ) +6 .D«-00*U ( 2 )-4 .D*-00*U ( 3)> 

1U(4) + (X=fr*4)*( A14-B( 1 ) 4-AA1*U(2)+AA2=!=U(2)-4=2) 

2 -q#x-*4 

RETURN 

15 0IFF8= U( 1 )-4.Df00X=U(2 )^■7.D+00*U(3)-4.D + 00*U(4)^■ 

1 ( X^’i'A) * ( Al«B ( 1 ) + AA 1«U ( 3 ) + AA2*U ( 3 ) =*^*2) 

2 -Q=<‘X=i'*4 
RETURN 

20 DIFF8= 2.D<-00*U( 2)-8.D«-00*U(3)+6.D+00*U(4M- 

1 (X«*4) *(A1>:^B(1)+AA1*U(4)+AA2=^U(4)**2) 

2 -q«x=?«4 
RETURN 

END 
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